Passive  Gust  Alleviation  for  a  Flying  Wing  Aircraft 


Dr.  Shijun  Guo 
Prof.  Otto  Sensburg 


Cranfield  University 
College  Rd 

Cranfield,  Bedford  MK43  OAL 
United  Kingdom 


EOARD  Grant  1 1-3073 

Report  Date:  January  2013 

Final  Report  from  26  August  2011  to  24  November  2012 


Distribution  Statement  A:  Approved  for  public  release  distribution  is  unlimited. 


Air  Force  Research  Laboratory 
Air  Force  Office  of  Scientific  Research 
European  Office  of  Aerospace  Research  and  Development 
Unit  4515  Box  14,  APO  AE  09421 


REPORT  DOCUMENTATION  PAGE 


Form  Approved  OMB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  the  burden,  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson 
Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply 
with  a  collection  of  information  if  it  does  not  display  a  currently  valid  OMB  control  number. 

PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


1.  REPORT  DATE  (DD-MM- YYYY)  2.  REPORT  TYPE  3.  DATES  COVERED  (From  -  To) 

10  January  2013  Final  Report  26  August  201 1  -  24  November  2012 


4.  TITLE  AND  SUBTITLE  5a.  CONTRACT  NUMBER 


Passive  Gust  Alleviation  for  a  Flying  Wing  Aircraft 


FA8655-1 1-1 -3073 

5b.  GRANT  NUMBER 


Grant  11-3073 

5c.  PROGRAM  ELEMENT  NUMBER 


6.  AUTHOR(S) 


Dr.  Shijun  Guo 
Prof.  Otto  Sensburg 


61102F 

5d.  PROJECT  NUMBER 

5d.  TASK  NUMBER 


5e.  WORK  UNIT  NUMBER 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Cranfield  University 
College  Rd 

Cranfield,  Bedford  MK43  0AL 
United  Kingdom 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

EOARD 

Unit  4515  BOX  14 
APO  AE  09421 


10.  SPONSOR/MONITOR’S  ACRONYM(S) 

AFRL/AFOSR/IOE  (EOARD) 


11.  SPONSOR/MONITOR’S  REPORT  NUMBER(S) 

AFRL-AFOSR-UK-TR-201 3-0008 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


14.  ABSTRACT 

This  final  report  presents  the  work  and  results  of  a  research  project  ‘Passive  Gust  Alleviation  for  a  Flying  Wing  Aircraft’  funded 
by  EOARD/US  AFRL  (Contract  FA8655-1 1-1-8073)  from  26  Aug.  2011  to  24  Nov.  2012.  In  the  project,  an  investigation  was 
made  into  the  technology  potential  of  a  passive  gust  alleviation  device  (PGAD)  and  its  application  to  a  Sensorcraft  of  high 
aspect  ratio  in  flying  wing  configuration.  It  is  aimed  at  minimizing  the  gust  response  of  the  aircraft  by  using  the  PGAD  integrated 
at  the  wing  tip.  The  project  has  been  carried  out  in  four  stages:  the  loading  analysis  including  aerodynamic  calculation  and 
mass  estimation,  structural  design  and  modeling,  gust  response  analysis  and  optimal  design  of  the  PGAD  for  minimum  gust 
response. 


Standard  Form  298  (Rev.  8/98) 
Prescribed  by  ANSI  Std.  Z39-18 


OARD/AFRL 


Cranfield 

I  UNIVERSITY 


2013 


Final  Technical  Report: 
(FA8655-1 1-1-3073) 


Passive  Gust  Alleviation  for  a  Flying  Wing  Aircraft 


Dr  Shijun  Guo 
Prof.  Otto  Sensburg 


Aerospace  Engineering 
Cranfield  University,  UK 


10.  Jan.  2013 


Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


UNIVERSITY 


2013 


EXECUTIVE  SUMMARY 


This  final  report  presents  the  work  and  results  of  a  research  project  ‘Passive  Gust  Alleviation 
for  a  Flying  Wing  Aircraft’  funded  by  EOARD/US  AFRL  (Contract  FA8655-1 1-1-8073)  from 
26  Aug.  2011  to  24  Nov.  2012.  In  the  project,  an  investigation  was  made  into  the  technology 
potential  of  a  passive  gust  alleviation  device  (PGAD)  and  its  application  to  a  Sensorcraft  of  high 
aspect  ratio  in  flying  wing  configuration.  It  is  aimed  at  minimizing  the  gust  response  of  the 
aircraft  by  using  the  PGAD  integrated  at  the  wing  tip.  The  project  has  been  carried  out  in  four 
stages:  the  loading  analysis  including  aerodynamic  calculation  and  mass  estimation,  structural 
design  and  modeling,  gust  response  analysis  and  optimal  design  of  the  PGAD  for  minimum  gust 
response. 

In  the  aerodynamic  analysis,  different  methods  including  vortex  lattice  and  CFD  methods 
were  employed  and  used  to  calculate  the  loading  at  specified  flight  cases.  In  the  preliminary 
design  phase,  mass  distribution  was  estimated  based  on  the  MTOW  of  55  tones  and  primary 
systems  including  the  structures,  airframe  systems,  power  plant,  fuel,  avionics  and  landing  gears. 
From  these  data,  the  shear  force,  bending  moment  and  torque  diagrams  acting  on  the  structure 
were  calculated. 

In  the  initial  layout  of  the  structure,  a  conventional  configuration  was  selected  with 
consideration  of  the  PGAD  at  the  wing  tip.  The  skin  panels  are  reinforced  by  I-section  stringers 
to  resist  buckling.  The  wing  is  mainly  made  of  carbon/epoxy  composites.  Structural  design 
including  stiffened  panel  buckling  and  stress  analysis  was  carried  out  by  using  the  FE  method  to 
make  sure  the  large  sweep  and  high  aspect  wing  structure  meets  the  design  requirements.  The 
initial  sizing  of  the  major  components  was  based  on  the  strength  and  buckling  criteria.  The  initial 
design  of  the  aircraft  structure  was  followed  by  detailed  stress  analysis  by  using  the  NASTRAN 
finite  element  method.  The  results  show  that  the  failure  indices  are  below  one  and  the  strains 
below  3600  ps  under  limit  load.  The  elastic  deflection  of  the  wing  at  semi-span  31.6  m  wing  tip 
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reaches  2.4  m  under  limit  load  in  the  worst  case.  The  flutter  speed  of  241  m/s  and  frequency  6.2 


Hz  for  the  full  fuel  case  was  predicted  for  the  aircraft. 

With  a  maximum  weight  and  a  span  of  63  m,  the  sensitivity  of  the  wing  to  the  gust  is  a  major 
concern.  Subsequently  modal  analysis  and  gust  response  of  the  wing  to  a  discrete  (1-cos)  gust 
load  in  a  range  of  equivalent  frequencies  were  calculated.  The  gust  model  is  in  compliance  with 
the  EASA  CS-25  airworthiness  requirement.  Two  flight  cases  of  zero-fuel  and  full-fuel  mass  at 
sea  level  were  considered.  Without  the  PGAD,  the  gust  response  amplitude  at  wing  tip  reaches 
approximately  4  m  for  the  worst  full-fuel  case  at  sea  level. 

Finally  the  gust  response  was  evaluated  by  taking  the  PGAD  with  optimized  key  design 
parameters  subject  to  normal  flight  constraints.  The  optimized  PGAD  together  with  the  wing 
aeroelastic  effects  is  very  effective.  For  the  case  where  rigid-body  motion  is  constrained  and 
PGAD  twist  angle  is  limited  to  10  degree,  the  gust  response  in  terms  of  wing  tip  deflection  and 
wing  root  bending  moment  can  be  reduced  by  over  18%  and  15%  respectively.  When  the 
transverse  rigid-body  motion  is  set  free  with  the  PGAD  twist  angle  limited  to  10  degree,  the  gust 
response  in  terms  of  wing  elastic  deflection  and  bending  moment  can  be  reduced  by  over  20% 
and  17%  respectively. 
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1.  Introduction 

This  project  is  aimed  at  minimizing  the  gust  response  of  a  flying  wing  aircraft  by  using  a 
passive  gust  alleviation  device  (PGAD)  mounted  at  the  wing  tip.  Although  various  passive  gust 
alleviation  technology  has  been  proposed  and  investigated  in  previous  research  [1-2]  especially 
for  similar  type  of  sensorcraft  [3-10],  this  proposed  PGAD  concept  presents  an  effective  option. 
The  concept  of  PGAD  is  illustrated  in  Fig. la  where  a  separate  rigid  wing  section  called  PGAD  is 
mounted  to  the  wing  tip  through  an  elastic  hinge.  The  elastic  hinge  is  made  of  a  torque  spring  and 
a  rotation  shaft,  which  is  mounted  to  the  front  spar  and  end  rib  of  the  wing  and  the  PGAD  as 
illustrated  in  Fig. la.  By  setting  the  rotation  shaft  axis  in  front  of  the  aerodynamic  center,  the 
PGAD  will  twist  nose  down  in  response  to  a  gust  load  to  alleviate  the  aerodynamic  force.  As 
illustrated  in  Fig.  lb,  the  gust  induced  wing  load  distribution  over  the  wing  would  be  significantly 
reduced  by  employing  the  PGAD.  For  this  particular  aircraft  of  a  large  sweep  back  wing,  the  wing 
bending-torsion  coupling  and  aeroelastic  effect  will  also  contribute  to  the  gust  alleviation.  As  the 
results,  it  can  minimize  the  negative  impact  of  gust  on  the  airframe  and  flight  performance. 


Figure  1.  (a)  PGAD  at  wing  tip,  (b)  Lift  distribution  with  and  without  PGAD 

The  nose  down  rotation  of  the  PGAD  surface  is  a  self-induced  passive  motion  in  response  to  an 
increase  of  aerodynamic  force  such  as  gust  load.  The  effectiveness  of  the  PGAD  depends  upon  the 
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key  design  parameters  and  the  wing  dynamic  behavior.  The  composite  wing  can  be  tailored  by 
optimizing  the  wing  structure.  The  primary  PGAD  design  parameters  are  the  torque  spring 
stiffness,  position  of  the  shaft  attachment  and  dimension  of  the  device.  The  torque  spring  stiffness 
and  position  of  the  shaft  attachment  determine  the  twist  angle,  quantity  of  gust  alleviation  and 
how  fast  the  PGAD  reacts.  The  dimension  of  the  device  scales  the  amount  of  gust  response 
alleviation  but  is  limited  by  the  flight  performance  in  normal  load  condition.  The  surface  area  on 
the  device  in  front  of  the  axis  "feels"  the  gust  first  and  tries  to  move  the  device  nose-up  in  a  very 
short  period,  but  soon  the  whole  device  will  then  rotate  nose  down  with  the  load  alleviated. 

To  develop  and  evaluate  the  PGAD  technology,  a  research  proposal  was  made  and  a  15-month 
contract  was  granted  by  EOARD/AFRL.  In  the  proposal,  a  work  programme  and  time  schedule  as 
shown  in  Table  1  was  set  for  the  project. 


Table  1.  Work  Programme  and  Time  Schedule 


Work 

Package 

Activity 

Elapsed  Time 
[months] 

i. 

Original  flying  wing 

i.i. 

Data  collection  and  FE  modelling  of  the  original  wing  structure  using  Patran  / 
Nastran  package 

2 

1.2. 

Fully  stressed  wing  structure  analysis  under  2.5 g  limit  load  and  3.75g  ultimate  load 
for  a  valid  wing  box  design. 

1 

1.3. 

Vibration,  flutter  and  divergence  analysis  of  the  baseline  design 

1.5 

1.4. 

Tuned  gust  stress  analysis  for  initial  evaluation  of  the  baseline  design. 

1 

2. 

Wing  model  with  gust  alleviation  device. 

2.1. 

Initial  analysis  to  determine  an  optimal  dimension,  torsion  spring  stiffness  and 
location  for  a  tuned  gust  under  the  design  constraint 

1.5 

2.2. 

Analysis  of  the  whole  wing  with  an  optimal  passively  rotating  section  for  maximum 
gust  alleviation  including  aeroelastic  effect 

1.5 

2.3. 

FE  model  of  the  wing  structure  with  an  integrated  gust  alleviation  device  -  an 
optimal  design  of  a  passively  rotating  wing  tip  section  for  a  tuned  gust 

1.5 

2.4. 

Vibration,  flutter  and  divergence  analysis 

1 

2.5. 

Tuned  gust  stress  analysis  at  wing  root  to  evaluate  the  device 

1 

3. 

Wing  model  with  soft  wing  tip. 

3.1. 

Determine  optimal  thickness  and  material  (laminate  layup)  for  the  tip  section 

0.5 

3.2. 

Tuned  gust  stress  analysis. 

0.5 

3.3. 

Vibration,  flutter  and  divergence  analysis. 

1 

4. 

Final  technical  report. 

1 

Total: 

15  months 
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A  large  aircraft  of  MTOW  55  tones  with  a  high  aspect  ratio  wing  span  63  m  and  flying  wing 
configuration  is  taken  in  the  case  study.  The  project  has  been  carried  out  in  four  stages  started 
from  the  loading  analysis  including  aerodynamic  calculation  and  mass  estimation  and  followed  by 
structure  design  and  analysis,  dynamic  and  gust  response  analysis  and  optimal  design  of  the 
PGAD  for  minimum  gust  response. 

In  the  Report  Chapter  2,  the  aerodynamic  analysis  results  by  vortex  lattice  and  CFD  methods  in 
specified  flight  cases  were  presented.  The  aircraft  mass  distribution  was  estimated  based  on  the 
MTOW  and  primary  systems  including  the  structures,  airframe  systems,  power  plant,  fuel, 
avionics  and  landing  gears.  From  these  data,  the  shear  force,  bending  moment  and  torque 
diagrams  acting  on  the  structure  were  calculated  and  presented. 

In  Chapter  3,  the  theoretical  base  for  analytical  analysis  of  the  aircraft  structure  was  presented 
first.  The  structure  initial  layout  in  a  conventional  configuration  was  selected  with  the  PGAD  at 
the  wing  tip.  Structural  design  including  stiffened  panel  buckling  and  stress  analysis  was  carried 
out  based  on  the  theory.  The  initial  design  of  the  aircraft  structure  was  followed  by  detailed  stress 
analysis  by  using  NASTRAN  finite  element  method  to  make  sure  the  structure  meet  the  design 
requirements.  The  results  show  that  the  failure  indices  are  below  one  and  the  strains  below  3600 
pe  under  ultimate  load.  The  wing  tip  deflection  of  the  wing  of  semi-span  31.6m  reaches  2.4  m 
under  limit  load  in  the  worst  case. 

In  Chapter  4,  modal  analysis  to  assess  the  structure  dynamic  behavior  was  carried  out  followed 
by  gust  response  analysis  of  the  wing  to  a  discrete  (1-cos)  gust  load  in  a  range  of  equivalent 
frequency  were  calculated.  The  gust  model  is  in  compliance  with  the  EASA  CS-25  airworthiness. 
Two  flight  cases  in  zero-fuel  mass  and  full-fuel  mass  at  sea  level  were  considered.  Without  the 
PGAD,  the  gust  response  amplitude  at  wing  tip  reaches  approximately  4  m  for  the  full-fuel  case  at 
sea  level. 

Finally  the  gust  response  was  evaluated  by  taking  the  PGAD  with  optimized  key  design 
parameters  subject  to  normal  flight  constraint.  The  optimized  PGAD  together  with  the  wing 
aeroelastic  effect  is  very  effective.  When  rigid-body  motion  is  constrained,  the  gust  response  in 
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terms  of  pure  elastic  deflection  of  the  wing  and  maximum  bending  moment  can  be  reduced  by 
over  18%  and  15%  respectively  with  the  PGAD  twist  angle  limited  to  10  degree.  When  the 
transverse  rigid-body  motion  is  set  free  with  the  same  PGAD  twist  angle  limit,  the  gust  response 
in  terms  of  wing  elastic  deflection  and  bending  moment  have  been  reduced  by  over  20%  and  17% 
respectively. 


2.  Aircraft  Load  Analysis 

2.1  Aircraft  data  and  aerodynamic  load 

The  aircraft  basic  data  given  is  listed  in  Table  2.  The  geometry  data  is  shown  in  Figure  A1  in 
Appendix  A.  The  cruise  speed  is  M=0.65  at  altitude  18.3  km.  The  flight  speed  for  critical  gust 
load  is  M=0.3  at  sea  level.  Without  the  detailed  wing  airfoil  data,  a  standard  NACA4415  is  chosen 
because  of  its  best  match  with  the  available  geometry. 


Table  2.  Design  technical  data  of  the  aircraft 


Wing 
semi 
span  (m) 

Fuselage 
length  (m) 

MTOM  (full 
fuel,  kg) 

MTOM  (empty 
fuel,  kg) 

Sweep 
angle  (deg) 

Cruise 
altitude,  km 

31.6 

14.7 

55350 

27674 

30 

18.3 

The  aerodynamic  pressure  and  load  distribution  over  the  3D  whole  aircraft  was  calculated  by 
using  CFD  method.  The  flow  velocity  and  vortex  plot  from  the  CFD  simulation  as  shown  in  Fig.3 
shows  the  source  of  the  pressure  drop  in  the  aircraft  body-wing  kink  region  and  the  wing  tip.  The 
aerodynamic  force  at  different  angle  of  attack  (AoA)  was  calculated  by  CFD  method.  Through  the 
analysis,  the  results  show  that  the  lifting  force  at  AoA=5  degree  meets  the  lift  requirement  in  full 
fuel  MTOM  case  at  cruise  M=0.65.  The  method  is  also  compared  with  vortex  lattice  method, 
panel  method  and  lifting  line  theory  with  the  spanwise  lifting  force  distribution  results  shown  in 
Fig.4.  From  the  results,  it  is  noted  that  the  3D  effect  of  wing  on  the  lift  distribution  along  the  span 
especially  the  outer  wing  captured  by  the  CFD  method  is  more  significant  than  the  other  methods. 
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Figure  3.  Flow  velocity  and  vortex  effect  on  the  pressure  at  the  wing  kink  region 


Figure  4.  (a)  Spanwise  lifting  force  and  (b)  aerodynamic  centre  by  different  methods 


2.2  Aircraft  mass  estimation 

The  aircraft  mass  distribution  was  estimated  based  on  the  MTOW  and  primary  systems 
including  the  structures,  airframe  systems,  power  plant,  fuel,  avionics  and  landing  gears.  Three 
fuel  tanks  are  located  in  the  wing  box  as  shown  in  Fig.5.  The  mass  and  location  of  the  engine  and 
landing  gears  are  also  shown  in  Fig.5.  The  system  mass  is  also  considered  in  the  example. 

Mass  distribution  for  half  aircraft  structures  is  5700  kg,  the  power  plant  and  landing  gears  are 
2500  kg,  systems  3500  kg  and  full  fuel  15900  kg  in  the  loading  evaluation.  Detailed  mass 
estimation  is  presented  in  Appendix  B.  The  resulting  shear  force  and  bending  moment  distribution 
along  the  span  is  calculated  and  shown  in  Fig.6. 
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Figure  5.  Half  aircraft  with  location  of  fuel  tanks,  engine  and  LG 


Figure  6.  (a)  Shear  force  and  (b)  bending  moment  diagram  in  different  cases 


3.  Analytical  and  Numerical  Methods 

3.1  Theoretical  study 
3.1.1  The  theory  and  analytical  method 

Since  no  structure  details  were  available,  an  initial  structure  analysis  was  carried  out  based  on 
the  calculated  load  shown  in  Chapter  2.  In  the  initial  design  stage,  the  structural  model  was 
simplified  to  a  thin-walled  composite  beam  model.  The  method  developed  by  Armanios  and  Badir 
[14]  and  the  dynamic  stiffness  method  [15]  were  used  and  described  below.  In  the  stiffness 
modeling  of  a  wing  box,  the  wing  was  divided  into  20  spanwise  segments  with  each  one  modeled 
as  a  uniform  thin-walled  double-cell  box  beam  between  the  leading  edge  and  rear  spar.  The  whole 
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wing  structure  was  modeled  as  an  assembly  of  those  box  beams  along  the  span.  A  relationship 


between  the  bending  moment  Mx,  torque  My  and  the  transverse  and  twist  deflections  at  the  end  of 
an  anisotropic  thin-walled  closed-section  beam  are  expressed  below. 


My  —  Cj2^P  +(^23^  and  M x  —  C23 <f)'  C y^h" 


(1) 


The  stiffness  coefficients  Cij  of  each  segment  can  be  calculated  based  on  its  geometry  and 
material  properties  and  integration  along  its  cross  sectional  circumference, 


(2) 


j 


where  Ae  is  the  enclosed  area  of  the  cross  section;  A(s),B(s)  and  C(s)  are  given  below. 


(3) 


In  the  above  equations,  A,y  is  the  coefficients  of  stiffness  matrix  (A)  of  the  composite  skin  and  spar 
webs  of  the  closed-section  beam.  According  to  the  force-deflection  relationships  in  Eq.  1  and 
stiffness  definition,  the  stiffness  coefficients  C33,  C22  and  C23  actually  represent  the  bending, 
torsion  and  bending-torsion  coupling  rigidities  of  the  wing  box  beam,  which  are  usually  expressed 
by  symbols  El,  GJ  and  CK  respectively.  Contribution  of  the  six  stringers  to  the  wing  box  bending 
stiffness  is  also  included  in  the  model. 

The  dynamic  stiffness  matrix  method  [15]  was  subsequently  used  for  the  dynamic  analysis.  In 
this  method,  the  equations  of  motion  for  each  of  the  thin-walled  box  beams  were  represented  as 
follows,  where  the  bending-torsion  stiffness  coupling  was  included  but  the  transverse  shear 
deformation  and  warping  effect  were  neglected. 


El  ■  h""  +  CK  •</>'"  +  m-h-m-  Xa -^  =  0 


(4) 
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GJ-f  +  CK-ti”  +  m-Xa-h-Ip<f>  =  0 


(5) 


where  h""  =  d4h/dy 4 ,  h  =  d2h/dt 2 ,  (j)m  =  d3h/dy3  and  (j  =  d2h/dt 2 .  By  solving  the  differential 


equations,  an  exact  solution  for  the  displacement  function  h(y)  and  <f>(y)  can  be  obtained.  A 


dynamic  stiffness  matrix  for  a  box  beam  can  be  subsequently  created  by  relating  the 
displacements  to  the  bending  moment  and  torque  at  both  ends  of  the  beam.  A  dynamic  stiffness 
matrix  for  the  whole  wing  box  structure  is  obtained  by  assembling  all  the  wing  box  beam  stiffness 
matrices  along  the  wing  span  direction. 

3.1.2  Gust  response  analysis 

It  is  noted  that  the  dynamic  stiffness  matrix  is  actually  a  combination  of  stiffness  and  mass 
matrices  of  the  beam  and  is  frequency  dependent.  Since  this  particular  type  of  matrix  produces  a 
non-standard  eigenvalue  problem,  it  is  solved  by  using  the  Wittrick-William  algorithm  [16],  By 
employing  the  normal  mode  method,  the  aeroelastic  equation  for  a  wing  coupled  with  shelf 
excited  unsteady  aerodynamic  forces  can  be  written  in  generalized  coordinates  as  follows.  The 
unsteady  aerodynamic  forces  were  calculated  by  using  the  classical  Theordorsen  theory  [17,18] 
and  the  strip  method  in  incompressible  airflow. 


[k„W] -\pV*[AL\,  +i4D\+i^pV-[AL\  {?}  =  0 


(6) 


With  gust  load  as  external  unsteady  aerodynamic  force,  the  aeroelastic  response  equation  of  the 
wing  structure  is  written  as 


[M]{x}  +  [D]{x}  +  [K]{x}  =  [ALl]{x}  +  [AL2]{x}  +  [ALi]{x}  +  {ALext} 


(7) 


where  [M],  [D],  [K]  are  the  structural  mass,  damping  and  stiffness  matrices;  [AL\\  and  [ALext]  are 
the  unsteady  aerodynamic  and  external  dynamic  force  matrices  respectively.  The  1  -cosine  discrete 
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gust  model  specified  in  the  airworthiness  regulation  is  used  in  this  current  investigation.  The  gust 
velocity  profile  of  a  1 -cosine  model  is  expressed  as  below. 


U  =  ^x[l-cos(f)]  (8) 

where  C/de  is  a  specified  design  gust  velocity,  s  is  the  distance  penetrated  in  the  gust,  and  H  is  gust 
gradient  distance  which  is  the  distance  parallel  to  the  airplane's  flight  path  for  the  gust  to  reach  its 
peak  velocity. 

3.1.3  Optimization  Method 

In  the  optimization  process,  the  gradient  based  determinant  method  (GBDM)  used  in  previous 
work  [19]  is  employed  to  determine  the  PGAD  design  variables.  Effort  is  primarily  focused  on 
minimizing  the  gust  response  and  loading  on  the  wing.  The  optimization  analysis  can  be 
expressed  as  follows: 


Minimise  F(x)  = 


R 


§  0 


\/\(Kv  K2,...dvd2...  )} 


(9) 

(10) 


where  f(x)  is  the  objective  function,  Rg(x)  wing  gust  response,  x  a  vector  containing  the  key 
parameters  of  the  PGAD  as  design  variables. 


3.2  Structural  layout  and  initial  analysis 

The  layout  of  the  primary  structure  of  the  flying  wing  aircraft  is  illustrated  in  Figure  7.  The 
structure  is  divided  into  11  spanwise  sections  in  the  modeling.  A  multi-spar  configuration  was 
chosen  for  the  inner  wing  and  a  conventional  two-spar  configuration  for  the  outer  wing.  For  the 
inner  wing,  the  front,  middle  and  rear  spars  are  located  at  14%,  50%,  and  80%  wing  chord.  For  the 
outer  wing,  the  front  spar  and  the  rear  spar  are  respectively  located  at  15%  and  75%  of  the  chord 
at  the  wing  tip  and  kept  parallel  to  the  leading  edge  along  the  span. 
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(a) 


Figure  7.  (a)  Structural  layout  and  (b)  sections  of  the  flying  wing  aircraft 


An  intermediate  modulus  carbon  fiber  epoxy  matrix  composite  (8552  epoxy  matrix  IM7  UD 
carbon  fiber)  has  been  chosen  for  the  wing  structure.  The  properties  of  the  material  are  presented 
in  the  Table  3.  Based  on  the  loading  and  structural  layout,  an  initial  structural  design  of  the  spars, 
ribs  and  skin  covers  was  carried  out.  The  results  of  the  initial  design  of  the  skin  panels  in  quasi¬ 
isotropic  layup  are  listed  in  Table  4. 


Table  3.  Some  technical  Data  for  aircraft  design  and  layout 


Ei 

(GPa) 

e2 

(GPa) 

G 

(GPa) 

V 

X, 

(MPa) 

xc 

(MPa) 

Y, 

(MPa) 

Yc 

(MPa) 

S 

(MPa) 

P 

(kg/m3) 

164 

12 

5.3 

0.32 

2724 

1690 

111 

246 

120 

1570 

Table  4.  Skin  panel  thickness  of  the  initial  design 


Section 

Upper  skin 
thickness  (mm) 

Lower  skin 
thickness  (mm) 

Section 

Upper  skin 
thickness  (mm) 

Lower  skin 
thickness 
(mm) 

i 

4.5 

3.7 

7 

5.2 

3.4 

2 

5.3 

2.9 

8 

4.7 

3.9 

3 

6.0 

3.1 

9 

4.2 

3.4 

4 

7.6 

5.2 

10 

3.1 

2.6 

5 

6.3 

4.5 

11 

2.1 

2.9 

6 

5.8 

3.9 

Prier  to  the  FE  modeling,  the  buckling  analysis  method  [11,12]  under  practical  design 
constraint  [13]  has  been  used  in  the  stiffened  skin  panel  design.  In  the  composite  structure 
analysis,  the  stress  level  was  limited  to  3500  micro  strain  under  damage  tolerance  constraint. 
Based  on  the  FE  model  of  the  structure,  which  satisfies  the  buckling  and  strength  requirements. 
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modal  analysis  and  initial  dynamic  response  to  a  discrete  gust  input  in  a  range  of  frequency  was 
conducted  without  considering  aeroelastic  effect. 

Based  on  the  initial  design,  a  buckling  analysis  for  the  stiffened  composite  skin  panels  is 
carried  out  [11-13].  As  an  example  result,  the  detailed  dimensions  of  the  upper  stringer-skin  panel 
in  the  kink  section  5  where  the  maximum  compression  load  occurs  are  given  in  Fig.  8.  The 
calculated  stringer  pitch  is  230  mm.  The  buckling  load  factors  of  wing  box  upper  surface  panels 
are  checked  with  ESDU  03001  [13].  The  overall  and  local  buckling  load  factors  for  this  panel  are 
2.1  and  1.1  respectively. 


K  ■ 

It 


£ 


T~ 


=F 


-  sm_STEINGER  FABEL  SIZING  results  - 

Stringer  Ease  Flange  ba  =  35.  21  mm  ta  =  4.  40  mm 

Stringer  Web  bw  =  68.40  mm  tw  =  3.42  mm 

Stringer  Top  Flange  bf  =  27.36  mm  tf  =  3.42  mm 

Stringer  Cross- section  Area  Ast  =  717.68  mm "2 

Skin  Cross-section  Area  Ask  =  1435.36  mm"2 

Stringer  Spacing  Fitch  =  228.27  mm 


Figure  8.  Detailed  dimensions  of  the  stringer-skin  panel  at  the  critical  section  5 


3.3  Structural  FE  model  and  stress  analysis 

In  this  investigation,  detailed  structural  stress  and  dynamic  analysis  of  the  wing  structure  made 

of  spars,  ribs  and  stringer  reinforced  skins  was  carried  out  by  using  the  FE  method  based 
NASTRAN  package.  In  the  FE  model,  the  skin,  ribs  and  spar  webs  are  modeled  by  using  shell 
elements;  the  stringers  and  spar  capes  are  modeled  as  beam  element.  In  the  dynamic  and  gust 
response  analysis,  the  engine,  LG,  fuel  tanks  and  control  devices  and  systems  were  modeled  as 
concentrated  mass  located  at  the  center  of  gravity  of  the  components. 

The  FE  analysis  results  show  a  maximum  displacement  of  2.56  m  at  the  wing  tip  under  limit 
load  (2.5g).  The  associated  failure  index  (FI)  of  the  upper  skin  is  under  0.51  as  shown  in  Fig.  9(a). 
The  maximum  strain  magnitude  reaches  3320  ps  which  is  under  the  limit  of  3600  ps  as  shown  in 
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Fig.  9(b).  It  can  be  observed  that  the  critical  strain  and  stress  are  located  around  the  kink  where 
the  stress  concentration  occurs  due  to  wing  geometric  change. 


Figure  9.  (a)  FE  results  of  FI  and  (b)  strain  plot  of  the  upper  wing  skin  under  limit  load 


The  structure  was  optimized  with  detailed  FE  analysis  presented  in  Appendix  E.  As  the  final 
results  under  limit  load,  the  maximum  strains  3570  ps  and  3560  ps  occur  in  the  lower  skin  and 
spars  respectively  (Appendix  E,  Fig.47).  The  results  indicate  that  the  preliminary  design  of  the 
wing  structure  meets  the  strength  requirements. 

3.4  Structural  FE  modal  and  gust  response  analysis 

Based  on  the  FE  model,  modal  analysis  was  carried  out  by  using  Nastran.  The  first  a  few 

frequencies  in  both  empty  fuel  and  full  fuel  cases  are  listed  in  Table  5,  with  the  associated  mode 
shapes  in  empty  fuel  case  presented  in  Fig.  10. 


Table  5.  The  first  a  few  frequencies  in  empty  fuel  and  full  fuel  cases 


Empty  mass 

Full-fuel  mass 

1st  bending  mode  (Hz) 

1.40 

0.69 

2nd  bending  mode  (Hz) 

5.96 

3.16 

3rd  bending  mode  (Hz) 

11.73 

6.80 

1st  torsion  mode  (Hz) 

16.04 

12.18 

Patran  2010.1  2  64-Blt  (MD  Enabled)  1 0-Aug-1 2 1 1 :30:27 

Fringe:  Default.  A1 6  Mode  1  :  Freq  =  1  3946.  Eigenvectors.  Translational.  Magnitude.  (NON-LAYERED) 
Deform  Default.  A16:Mode  1  Freq.  =  1 .3946.  Eigenvectors.  Translational. 


1.64+000 
1  44+000 
1  34+000 
1.24+000 
1.13+000 
1.03+000 


9.27-001 
8  24-001 
7.21-001 
6.18-001 
6.16-001 
4.12-001 
3  08001 
2.06-001 


Patran  201 01  2  64-Brt(MO  Enabled)  KMug-12 11  32  41 

Fringe  Default  A16  Mode  1 1  Freq  ■  6  9642  Eigenvectors  Translational  Magnitude  (NON-LAYERED) 
Deform  Default  A16  Mode  1 1  Freq  •  6  9642  Eigenvectors  Translational, 


1  66+000 
1  46+000 
1  36+000 
1  24+000 
1  14+000 

1  04+000 
9  32-001 
828601 
7  26601 
621-001 
618601 

4  14-OOtU 
3  1  160)1 

2  07-00)1 
1  04-00 1 B 
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Palran  2010 1  2  04-Btl  (MD  ErvabM  I0-Aug-12  1 1  33  31 

Fringe  OefaKt  A16  Mode  31  Freq  •  11  73  Eigenvectors,  Translational  Magnitude.  (NON-LAYERED) 
Deform  Detaiit  A16  Mode  31  Freq  •  1 1 .73.  Eigenvectors  Translational 


9  88-001 
9  23-001 
8  67-001 
791-001 
7  26-001 
6  69-001 
6  93-001 
6  27-001 
4  61-001 
3  96-001 
3  29-001 
264-001 
1  96-001 
1  32-00l| 
6  69-005 


r 

1 

I 

1 


I  08*000 
101*000 
9  34-001 
8  62-001 
7  90-001 


7  18-001 
6  47-001 
6  76-001 
6  03-001 
4  31-001 
3  69001 
2  87-001 
2  16-001 
I  44-001 1 


Figure  10.  Normal  modes  of  the  wing  in  empty  mass  case 


In  the  gust  response  analysis,  three  values  of  gust  gradient  were  selected  to  consider  the  whole 
range  from  9  m  to  107  m.  The  three  values  and  equivalent  gust  frequencies  are  shown  in  Table  6. 
One  of  them  is  the  typical  gust  gradient  of  12.5  times  the  mean  chord.  It  is  noted  that  the  gust  load 
is  in  the  range  of  the  first  three  bending  modes  of  the  wing  structure.  This  causes  a  concern  of  the 
wing  structure,  which  is  likely  to  be  sensitive  and  have  large  response  to  gust  load.  The  flight 
cases  and  data  considered  in  the  structural  and  gust  load  analysis  are  listed  in  Table  7. 


Table  6.  The  gust  gradient  values  and  equivalent  frequency 


H-9  m 

#=12.5*chord=79.1  m 

tf=107 m 

Frequency  (Hz) 

5.67 

0.65 

0.48 

Table  7.  Difference  flight  cases  for  gust  response  analysis 


Gust  Case  1 

Gust  Case  2 

Cruise  Case  1 

Cruise  Case  2 

Full  Fuel 

Empty 

Full  Fuel 

Empty 

Altitude  (ft) 

Sea  level 

Sea  level 

60000ft 

60000ft 

Weight  (kg) 

27674.22 

11729.5 

27674.22 

11729.5 

Mach  No 

0.3 

0.3 

0.65 

0.65 

Gust  Load  Factor 

2.95 

3.81 

2.13 

3.42 

In  the  initial  evaluation  of  the  gust  response  based  on  the  FE  model,  the  gust  load  was 
calculated  and  applied  as  an  external  dynamic  force  on  the  aircraft  clamped  at  the  body  center  line 
in  the  same  spanwise  distribution  as  the  aerodynamic  force.  Structural  damping  was  considered, 
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but  the  rigid  body  motion  is  constrained.  The  aerodynamic  damping  and  aeroelastic  coupling  was 
ignored.  The  most  critical  case  at  M=0.3  and  sea  level  was  considered  in  this  case  study. 

In  the  empty  fuel  case,  the  gust  response  measured  as  displacement  at  wing  tip  is  shown  in 
Fig.  11  (a).  From  the  results,  it  is  noted  that  the  maximum  response  was  due  to  the  gust  of  larger 
gust  gradient  H= 79m  and  107m  corresponding  to  lower  frequency  range.  Although  the  gust 
frequency  is  below  the  1st  bending  mode,  it  is  likely  to  excite  the  1st  mode  and  produce  large 
response.  It  is  also  noted  that  the  gust  frequency  of  5.67  Hz  corresponding  to  the  gust  gradient 
H=  9m  is  very  close  to  the  2nd  mode  of  5.96  Hz  and  causes  a  long  oscillation.  However  the 
deflection  is  much  smaller  than  the  larger  gust  gradient  cases. 

In  the  full  fuel  case,  the  gust  response  is  shown  in  Fig.  11(b).  From  the  results,  it  is  noted  that 
the  maximum  response  was  corresponding  to  the  gust  gradient  H= 79m  with  the  frequency  0.65  Hz 
very  close  to  the  wing  1st  bending  mode  of  0.69Hz.  This  is  the  most  critical  case  to  be  further 
investigated  by  considering  the  PGAD. 


(a) 


(b) 


Figure  11.  The  gust  response  (a)  in  the  empty  fuel  case  and  (b)  in  the  full  fuel  case 


4.  Gust  Response  Analysis  of  the  Wing  with  PGAD 

4.1  Beam  model  analysis 

4.1.1.  Wing  beam  model  and  gust  response 

The  wing  structure  is  then  simplified  by  using  beam  model  to  consider  the  unsteady 
aerodynamics  and  aeroelastic  coupling  in  the  gust  response  analysis  based  on  the  theory  presented 

14 

Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


2013 


EOARD/AFRLl 


Cranfield 

I  UNIVERSITY 


in  section  2.  The  beam  model  is  divided  into  20  sections  along  the  wing  span  plus  one  section 
representing  the  PGAD  device  of  1.85m  in  length  connected  to  the  beam  model  by  a  rotational 
spring.  The  nodes  are  located  at  shear  center  of  wing  sections  assumed  at  40%  of  the  chord.  The 
section  bending  and  rotational  stiffness  are  taken  from  cross  sections  of  the  section  planes 
perpendicular  to  the  neutral  axis  defined  by  21  nodes.  The  mass  distribution  is  transferred  from 
the  FE  model.  This  makes  sure  that  the  same  bending  and  rotational  stiffness  and  similar  mass 
distribution  as  the  FE  model  are  used  for  the  beam  model.  Table  8  shows  that  the  dynamic 
behavior  of  the  simplified  beam  model  in  the  full  fuel  case  is  very  close  to  the  FE  model. 


Table  8.  The  modal  frequency  from  the  wing  FE  and  beam  models 


FE  Model 

Beam  Model 

1st  Bending  Mode: 

0.69  Hz 

0.626  Hz 

2nd  Bending  Mode: 

3.16  Hz 

3.104  Hz 

3rd  Bending  Mode: 

6.80  Hz 

7.500  Hz 

1st  Torsion  Mode: 

12.18  Hz 

13.22  Hz 

In  the  following  gust  response  analysis,  only  the  most  critical  full  fuel  at  sea  level  case  was 
considered.  Fig.  12(a)  shows  the  gust  response  of  the  wing  to  the  three  gust  load  frequencies  listed 
in  Table  6  without  the  PGAD.  Fig.  12(a)  also  shows  the  maximum  response  occurring  in  the 
typical  gradient  length  at  0.65  Hz.  Then  the  PGAD  with  three  different  design  parameters  was 
considered  and  results  are  shown  in  Fig.  12(b)  in  comparison  with  the  response  without  PGAD.  In 
the  figure,  a=-0.3  and  -0.9  represents  the  location  of  the  rotation  shaft  is  mounted  at  30%  and  90% 
of  the  semi-chord  in  front  of  the  mid  chord  of  the  PGAD.  The  torque  spring  stiffness  considered 
in  the  initial  study  is  40  kNm  and  10  kNm. 

Figure  12(b)  shows  that  the  gust  response  for  the  device  shaft  location  a=-0.3  (located  at  35% 
chord  from  leading  edge)  is  greater  than  the  case  without  the  device.  This  is  due  to  the  shaft 
location  being  behind  the  pressure  center  of  the  surface  and  causing  a  positive  rotation  and  lift. 
While  the  gust  response  for  the  device  shaft  location  a=-0.9  (located  5%  chord  from  leading  edge) 
is  significantly  reduced  especially  for  the  lOkNm  softer  spring  case.  The  results  in  Fig.  12(b) 
shows  that  for  the  case  of  a=  -0.9  and  torque  spring  stiffness  40kNm,  the  gust  response  decrease 
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by  30%  from  2.2m  to  1.5m.  This  is  due  to  a  negative  pitching  moment  generated  by  the  device 
when  the  neutral  axis  is  located  ahead  of  the  aerodynamic  center. 


(a)  (b) 

Figure  12.  The  gust  response  of  the  wing  (a)  without  PGAD  and  (b)  with  the  PGAD 


The  device  was  set  in  a  pitching  degree  of  freedom  from  +10  degree  to  -20  degree  at  the 
connection  section.  In  response  to  the  gust,  the  wing  deflects  and  rotates  along  the  span  and 
produce  unsteady  aerodynamic  forces.  Fig.  13(a)  shows  the  wing  angle  of  attack  (twist  angle) 
associated  with  the  gust  response  shown  in  Fig.  12(b)  at  the  section  where  the  PGAD  is  mounted. 
Fig.  13(b)  shows  the  sum  of  the  PGAD  twist  angle  and  the  wing  AoA  at  the  wing  tip.  The  gust 
generates  a  negative  structure  twist  angle  as  a  consequence  of  bending  and  torsion  coupling  of  the 
large  sweep  back  wing. 


(a) 


(b) 


Figure  13.  (a)  AoA  of  the  wing  at  wing  tip  (b)  twist  angle  of  the  PGAD 
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For  a  low  rotational  stiffness  spring,  the  gust  generates  a  negative  twist  angle  on  the  device  and 
causes  a  sudden  load  reduction  with  a  negative  displacement  at  the  early  stage.  With  the  increase 
of  gust  induced  angle  of  attack,  the  lift  and  the  response  increase.  The  low  spring  stiffness  allows 
the  device  twist  reaches  -20  degree  limit  as  displayed  in  Fig.  13(b). 

By  comparing  the  amplitude  of  different  neutral  axis  location  and  spring  rotational  stiffness 
during  the  whole  gust  process,  an  optimal  shaft  location  and  spring  stiffness  can  significantly 
reduce  the  gust  response. 

4.1.2.  PGAD  optimization  for  minimum  gust  response 

The  objective  for  optimization  is  to  minimize  the  gust  response  by  varying  the  shaft  location 

and  torque  spring  stiffness.  A  gradient  based  optimizer  is  applied  in  this  optimization  process. 


MinV^x') 


—  0.  1  >  a  >  —0.  9 
80k Nm>  Ks  >  10k A/m 


(11) 


The  torque  spring  stiffness  is  set  as  the  design  parameter  at  a  given  shaft  location.  The  spring 


stiffness  at  three  different  shaft  locations  a=  -0.1,  -0.7  and  -0.9  are  optimized  to  achieve  a 


minimum  response.  The  shaft  location  at  a=  -0.5  which  is  the  aerodynamics  center  location  is 
taken  as  reference. 

Fig.  14(a)  shows  the  optimal  spring  stiffness  for  the  three  shaft  locations  are  80kNm,  15kNm, 
and  31kNm  with  gust  response  of  2.41m,  1.11m,  1.2m  respectively.  Comparing  with  the  reference 
response,  the  optimized  stiffness  leads  to  the  gust  response  reduced  by  48%  and  45%  for  a=-0.7 
and  -0.9.  The  response  decreases  along  with  the  shaft  moving  forward  and  spring  stiffness 
decreasing. 

A  slightly  smaller  response  occurs  at  a=  -0.7  rather  than  a  =  -0.9.  This  is  due  to  that  a  softer 
spring  is  quicker  and  easier  to  reach  the  lower  bound  of  the  twist  angle  limit  for  the  device. 
However  a  sudden  lift  reduction  could  result  in  a  higher  response  due  to  aeroelastic  effect  from 
the  wing  structure  bending-torsion  coupling. 
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(a)  (b) 

Figure  14.  (a)  Optimized  response  varying  with  the  spring  stiffness 
(b)  Optimization  history  for  different  shaft  locations 


In  order  to  maintain  the  normal  flight  performance  without  gust,  a  practical  design  constraint  is 
set  to  limit  the  PGAD  twist  angle  to  -2  degree  during  cruise.  This  will  prevent  a  significant  lift  lost 
in  normal  flight.  For  the  case  of  a=-0.7,  the  lower  bound  of  the  spring  stiffness  is  28kNm.  Under 
this  limit,  the  gust  response  is  1.79m  with  a  gust  load  reduced  by  17%  comparing  with  the 
reference  case. 

4.2  3D  FE  model  analysis 

Based  on  the  3-D  FE  model  of  the  whole  aircraft  (details  in  Appendix  E),  further  gust  analysis 
was  carried  out  with  the  PGAD  spanwise  dimension  increased  from  the  previous  1.6  m  to  2  m. 
However  the  torque  spring  stiffness  was  designed  to  limit  the  PGAD  twist  angle  within  10  degree. 
To  minimize  the  flow  turbulence  between  wing  and  PGAD  in  action,  the  interface  is  set 
streamwise  as  shown  in  Fig. la.  Accordingly  the  shaft  should  be  normal  to  the  interface  to  make 
sure  the  PGAD  can  rotate  freely.  However,  the  PGAD  could  be  designed  with  its  elastic  axis 
either  in  parallel  to  the  shaft  or  to  the  front  spar.  Both  cases  were  studied  with  the  1st  case  results 
presented  in  the  following  section  and  the  second  less  effective  case  presented  in  Appendix  F.  The 
shaft  location  in  chordwise  (a)  and  spring  stiffness  ( GJ)  were  taken  as  independent  design 
variables.  In  the  FE  modeling,  the  shaft  and  torque  spring  are  modeled  by  using  spring  elements. 
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The  following  gust  alleviation  analysis  is  presented  in  two  parts.  Firstly  the  rigid  body  motion 
of  the  aircraft  was  restricted  by  clamping  the  aircraft  body  center  line  as  previous  study.  Then  the 
aircraft  is  set  free-free  to  take  the  rigid  body  motion  effect  into  account. 

4.2.1.  Gust  response  for  the  shaft  in  parallel  to  global  Y-axis 

In  the  four  different  locations  of  the  shaft  a=-0.3,  -0.5,  -0.7  and  -0.9  studied  in  previous  beam 

model,  only  the  practical  and  effective  locations  a=-0.3  and  -0.7  with  the  typical  gust  gradient 
H=  12.5  times  chord  was  taken  in  the  analysis. 

In  the  case  a=-0.7,  the  shaft  is  located  exactly  on  the  front  spar  (15%  chord  from  leading  edge), 
which  is  a  practical  position.  A  series  of  spring  GJ  values  were  studied  and  3  typical  results  were 
presented  in  Table  9.  Comparing  with  the  previous  results,  much  greater  gust  alleviation  over  35% 
can  be  achieved  when  the  PGAD  twist  angle  was  limited  to  -20  degree.  Under  the  practical 
constraint  for  the  PGAD  twist  angle  -10°  set  as  boundary,  over  18%  gust  alleviation  and  15% 
bending  moment  (wing  root)  reduction  have  been  achieved.  The  torque  spring  stiffness  was 
increased  to  58  KNm/rad.  Taken  this  spring  stiffness,  the  gust  response  results  in  terms  of  PGAD 
twist  angle,  wing  tip  displacement  and  wing  root  bending  moment  for  different  shaft  locations  are 
shown  in  Fig.  15-17.  Some  of  the  results  have  been  published  in  a  conference  [31]. 


Table  9.  Gust  response  with  PGAD  shaft  normal  to  streamline 


Case 

Spring 

stiffness 

(Nm/rad) 

Wing  tip 
Disp. 

(m) 

Disp. 

Reduction 

Bending 

Moment 

(KNm) 

BM 

Reduction 

PGAD 

Relative 

Twist 

angle(°) 

Initial 

design 

/ 

3.07 

/ 

5600 

/ 

/ 

2.8E4 

1.99 

35.2% 

3987 

28.8% 

-19.7 

a=-0.7 

4.0E4 

2.28 

25.7% 

4461 

20.3% 

-14.2 

5.8E4 

2.51 

18.2% 

4760 

15% 

-10.0 

a=-0.5 

5.8E4 

2.59 

15.6% 

5240 

6.4% 

-8.3 
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Figure  15.  PGAD  twist  angle  in  response  to  gust  Figure  16.  Wing  tip  displacement  response 


Figure  17.  Wing  bending  moment  at  root 


All  the  above  analysis  produces  relative  values  of  the  gust  response  to  assess  the  effectiveness  of 
the  PGAD.  Since  the  aircraft  was  clamped  at  its  root,  all  the  input  energy  from  gust  is  assumed  to 
be  absorbed  by  the  elastic  structure.  This  results  in  an  overestimation  of  the  gust  response  and 
loading  to  the  structure.  In  order  to  predict  the  real  lift  gust  response,  the  constraint  set  to  the 
aircraft  should  be  removed  to  allow  for  a  free-body  motion  as  real  life.  The  analysis  and  results 
are  presented  in  the  following  section. 

4.2.2.  Gust  response  of  the  whole  aircraft  with  free-body  motion 

In  this  free-body  aircraft  case  study,  only  the  degree  of  freedom  in  Z-direction  was  set  free  so 

that  the  wing  can  plunge  freely  in  transverse  direction.  The  reason  for  not  being  able  to  remove 
the  pitching  constraint  was  because  the  whole  aircraft  was  not  trimmed  for  stable  level  flight  yet. 
The  aircraft  flight  stability  is  beyond  the  current  project  scope.  Nevertheless  the  transverse  free 
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rigid-body  motion  makes  the  analysis  condition  much  closer  to  real  lift.  In  the  analysis,  only  the 
most  practical  and  effective  shaft  location  a=-0.7  was  taken  as  example.  As  shown  in  Figure  18, 
the  PGAD  twist  angle  is  less  than  -10°  (-0.175  rad)  within  the  limit.  The  aircraft  body  (CG 
position)  displacement  and  the  wing  elastic  deformation  measured  at  tip  and  the  combined  elastic 
and  rigid  body  displacement  at  wing  tip  are  presented  in  Fig.  19.  To  compare  with  the  previous 
results,  the  wing  pure  elastic  gust  response  with  and  without  PGAD  is  shown  in  Fig.  20.  Figure  21 
shows  the  resulting  bending  moment  measured  at  wing  root  in  response  to  the  gust  load.  The 
results  show  that  the  gust  induced  maximum  wing  elastic  deflection  0.65m  without  PGAD  is 
reduced  by  20%  to  about  0.52m  and  the  bending  moment  is  reduced  by  17%. 


Figure  18.  PGAD  relative  twist  angle  Figure  19.  Wing  tip  displacement  comparison 
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Figure  21  Wing  bending  moment  response 
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5.  Conclusions 

The  PGAD  design  concept  was  evaluated  in  this  investigation.  A  full  scale  flying  wing  aircraft  of 

high  aspect  ratio  wing  has  been  taken  as  an  example  to  demonstrate  the  effectiveness  of  the 

PGAD  for  gust  alleviation.  Through  the  study,  the  following  conclusions  can  be  drawn. 

•  The  aerodynamic  analysis  has  been  carried  out  by  using  high  fidelity  modeling  methods. 
Detailed  pressure  distribution  and  flow  features  especially  in  the  kink  and  wing  tip  regions 
have  been  predicted  by  using  CFD  method  to  produce  accurate  loading  results. 

•  The  mass  distribution  of  the  whole  aircraft  with  the  primary  structure,  systems,  components 
and  payload  has  been  estimated  for  load  prediction.  However  the  mass  distribution  and  CG 
location  has  not  been  tuned  to  meet  level  flight  trim  condition  since  the  flight  stability  is 
beyond  the  project  scope. 

•  Based  on  an  initial  design,  detailed  stress  analysis  of  the  whole  aircraft  structure  modeled  by 
FE  method  based  on  NASTRAN  has  been  conducted  to  meet  the  design  requirements.  The 
maximum  strain  of  the  composite  structure  is  limited  to  3600  ps  considering  damage  tolerance. 

•  Following  the  modal  analysis  of  the  structure,  flutter  speed  of  241  m/s  has  been  predicted  using 
V-g  method  by  NASTRAN. 

•  Since  the  fundamental  frequency  of  this  particular  high  aspect  ratio  flexible  wing  structure  is 
only  0.69  Hz,  the  aircraft  is  very  sensitive  to  the  gust  in  the  whole  range  of  gust  length.  A 
significant  gust  response  occurs  in  the  specified  flight  condition  at  sea  level.  The  investigation 
shows  that  a  significant  reduction  of  gust  response  can  be  achieved  by  using  the  PGAD.  By 
optimizing  the  PGAD  key  parameters,  a  minimum  gust  response  can  be  achieved. 

•  When  the  rigid-body  motion  of  the  aircraft  is  constrained,  up  to  35%  gust  alleviation  can  be 
achieved  by  using  the  optimized  PGAD  subject  to  -20°  twist  angle  limit.  When  the  twist  angle 
is  limited  to  -10  °,  the  gust  response  in  terms  of  wing  tip  deflection  and  wing  root  bending 
moment  can  be  reduced  by  over  18%  and  15%  respectively. 
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•  When  the  transverse  rigid-body  motion  is  set  free  under  the  10°  twist  angle  limit,  the  gust 
response  in  terms  of  wing  elastic  deflection  and  bending  moment  can  be  reduced  by  over  20% 
and  17%  respectively. 
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Appendix  A.  Aerodynamic  Analysis 

A.l  Initial  analysis  and  Modified  wing  tip 
A.1.1  Wing  Geometry 

The  original  geometry  of  the  wing  was  presented  in  Fig.Al  as  following: 

All  coordinates  in  mm 


Figure  A1 .  Initial  geometry  of  the  wing 

In  order  to  carry  the  aerodynamic  study,  an  aerofoil  had  to  be  chosen  to  fit  as  close  as  possible  to 
coordinates  given  in  the  specifications.  The  only  data  given  on  the  aerofoil  shape  was  the 
thickness  of  15%  of  the  chord.  A  research  of  the  available  aerofoils  for  the  specified  thickness  was 
conducted.  From  the  sixteen  aerofoils  selected,  none  was  fitting  perfectly  with  the  coordinates  of 
the  geometry.  However,  with  the  approval  of  Dr  Guo,  the  standard  aerofoil  NACA  4415  as  shown 
in  Fig.  A.2  was  chosen,  as  it  has  one  of  the  highest  lift  coefficients  for  AoA=0°  (-0.50  for 
Reynolds  Numbers  of  106).  Indeed,  zero  angle  of  attack  was  the  reference  angle  of  attack  to  begin 
the  aerodynamic  analysis. 
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Figure  A2.  NACA  4415  aerofoil 


A.1.2  Flight  Cases 

Two  different  cases  were  used  within  the  aerodynamic  analyses  and  then  compared  [21]: 

•  Critical  Gust  Case  at  Full-Fuel  Weight,  Sea  level,  Mach  0.30: 

Air  density  p=  1.225  kg/m3 
Kinematic  density  v=  1.46x1 0-5  m2/s 
V elocity  v=  1 02  m/s 
Mass:  27674.22  kg  (semi  span) 

•  Cruise  Case  at  Full-Fuel  Weight,  60000  ft,  Mach  0.65: 

Air  density  p=  0.115  kg/m3 
Kinematic  density  v=  1.23x1 0-4  m2/s 
Velocity  v=191.8  m/s 
Mass:  27674.22  kg  (semi  span) 

The  aerodynamic  analyses  of  these  study  cases  were  first  performed  with  zero  angle  of  attack. 

A.2.1  Initial  CFD  analysis 

Initial  calculations  of  the  lift  distribution  of  this  geometry  were  conducted  by  using  the  CFD 
software  FLUENT.  From  the  first  results  of  the  analysis,  it  appeared  a  sudden  drop  of  the  lift  at 
the  tip  of  the  wing.  As  shown  in  the  Fig.  A3,  the  lift  decreases  significantly  and  becomes  negative 
from  30m  spanwise  to  the  wing  tip.  This  phenomenon  can  be  explained  by  the  turbulences  created 
by  the  angled  wing  tip.  Indeed,  the  wing  tip  shape  has  a  large  influence  on  the  aerodynamic 
characteristics  at  the  tip.  Investigations  have  been  conducted  on  the  influence  of  the  tip 
characteristics  on  the  wing  aerodynamic  characteristics  [20].  This  study  has  shown  by  wind  tunnel 
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experiments  that  the  non-alignment  of  the  tip  boundary  with  the  streamlined  flow  increases  local 
drag  and  lowers  the  local  lift. 

Spanwise  distribution  of  lift 

Wing  with  angled  tip,  final  results 

25000 


Spanwise  position  [m] 


Figure  A3.  Lift  distribution  for  the  geometry  with  angled  tip  (AoA=0°,  M=0.30) 

A.2.2  Modified  wing  tip  geometry 

The  original  angled  wing  tip  geometry  will  affect  the  accuracy  of  CFD  simulation.  The  loss  of 
pressure  would  also  have  negative  effect  on  the  PGAD  effectiveness.  Therefore,  a  modified 
geometry  was  proposed  to  minimize  the  impact.  The  original  wing  was  extended  and  then  cut  at 
y=31.14  m  in  parallel  to  xz  plane  to  make  the  tip  section  aligned  with  the  upstream  flow  velocity 
as  shown  in  Fig.  A.4 
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All  coordinates  in  mm 
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Figure  A4.  Final  wing  geometry 

With  this  new  configuration,  the  CFD  analysis  showed  that  the  drop  of  the  lift  at  the  tip  was 
eliminated.  With  this  tip  configuration,  the  lift  distribution  has  a  smoother  shape  at  the  tip  and 
keeps  a  positive  value  at  the  tip. 

A.2  Aerodynamic  methods 

A.2.1  The  Lanchester-Prandtl  Theory 

The  lift  distribution  of  a  three-dimensional  wing  is  not  the  successive  analyses  of  the  two- 
dimensional  cross-sections  along  the  span.  Indeed,  the  lift  created  on  each  wing  section  depends 
on  the  characteristics  of  the  section  but  also  on  the  lift  created  by  the  neighbouring  sections.  The 
Lanchester-Prandtl  Theory,  named  also  as  the  Prandtl’s  Lifting  Line  Theory  is  a  theoretical  model 
of  the  lift  distribution  based  on  the  three-dimensional  wing  geometry  [22],  This  theory  was 
elaborated  at  the  beginning  of  the  20th  century  and  is  still  used  for  preliminary  calculations. 

The  main  assumptions  of  this  theory  are: 

The  flow  is  a  steady  potential  flow:  inviscid  and  irrotational 
The  flow  is  incompressible 
High  aspect  ratio  wing 
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-  Low  sweep  angle 

The  three-dimensional  wing  is  discretised  into  multiple  span  wise  sections.  Each  section  is 
modelled  by  a  horseshoe  vortex,  which  is  the  combination  of  a  bound  vortex  and  two  free-trailing 


vortices  which  extend  downstream.  Multiple  bound  vortices  are  all  aligned  along  a  single  line,  the 


lifting  line  represented  in  the  Figure  A5,  conventionally  located  at  the  quarter  chord  line. 


Figure  A2.  Horseshoe  vortices  and  the  Lifting  Line  [22] 


The  theory  relates  the  circulation  r  and  the  lift  per  unit  span  using  the  Kutta-Joukowski  theorem: 

L'  =  Poo  T  (A.l) 

Where: 

L'  is  the  lift  per  unit  span 

Poo  is  the  freestream  density  of  the  fluid 

Vqo  is  the  freestream  velocity  of  the  fluid 

The  downwash  induced  by  the  trailing  edge  vortices  is  calculated  for  each  horseshoe  vortex.  The 
downwash  at  the  point  y0  is  induced  by  all  the  trailing  vortices  along  the  lifting  line.  The  trailing 
vortex  at  the  coordinate  y  creates  a  downwash  dw  at  the  point  y0  given  by  the  Biot-Savart  law: 


(dr /dy)dy 

4<yo  -  y) 


(A.2) 


The  total  downwash  at  the  point  y0  is  the  sum  of  the  effects  of  all  the  trailing  vortices  along  the 
lifting  line: 
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l  rb/2 


(dr /dy)dy 


(A.3) 


-6/2  yo-y 

The  downwash  affects  the  local  flow  along  the  span.  The  effective  angle  of  attack  cteff  depends 


on  the  geometrical  angle  of  attack  a  of  the  section  but  also  on  the  induced  angle  of  attack  due  to 
the  downwash  af. 

«e//(yo)  =  a(y0)  -  «i(yo)  (a.4) 

Where: 


«i(y0) 


1  rb/2  (d£/dy)dy 
4nVmJ_b/2  y0-y 


(A.5) 


And 


<*eff(yo) 


r(y0) 

nVmc(y0) 


+  <*l= o 


(A.6) 


aL=0  is  the  zero-lift  angle  of  attack,  known  from  the  aerofoil  characteristics 


Figure  3:  Effect  of  the  downwash  on  the  angle  of  attack  [22] 


The  unique  unknown  of  these  equations  is  the  circulation  r .  To  calculate  it,  the  following 
transformations  are  done: 


And 


b  ,  , 

y  =  —  —  cos  0  (A. 7) 

r(0)  =  ro  sin6>  (A. 8) 
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With  0  <  9  <  n 

The  Fourier  series  are  used  to  approximate  the  expression  of  the  circulation.  The  accuracy  of  the 
approximation  depends  on  the  number  of  terms  used.  The  more  terms  are  used,  the  more  accurate 

is  the  approximation,  but  the  more  complex  will  be  the  resolution  of  the  system. 

So, 


N 


r(0o)  =  2bVm  ^  An  sin  nd0  (A.9) 


Hence,  the  equation  (A.4)  becomes: 


2b  n  n 

a(0°)  =  nc^0  ^  ^  K  sin  n60  +  aL=0 (0O)  +  ^  n  An 


sinn0o 

sin0n 


(A.  10) 


i  i 

For  a  given  00,  which  corresponds  to  a  spanwise  location,  the  previous  equation  is  evaluated. 
Hence,  as  there  are  N  An  unknowns,  the  procedure  is  repeated  for  N  different  spanwise  locations 
to  be  able  to  solve  the  system. 

From  these  equations,  the  drag,  the  local  lift  coefficient  can  also  be  calculated. 

In  order  to  compare  the  different  methods  of  calculation  of  the  lift  distribution,  an  equivalent 
model  of  the  wing  studied  was  created.  A  straight  rectangular  wing  of  a  semi-span  of  31.14  m, 
chord  of  4.57  m,  using  NACA  4415  was  analysed.  This  wing  model  is  based  on  the  wing 
geometry  in  the  specifications.  The  outer  part  of  the  wing,  from  the  kink  at  10.5  m  to  the  tip  is 
taken  as  a  model.  The  swept  back  angle  is  neglected  in  this  application  and  the  same  geometry  as 

from  the  outer  part  is  used  for  the  inner  wing  geometry  from  the  root  to  the  kink. 

A 


Chord  =  4.57  m 


Aerofoil  NACA  4415,  AoA=0° 
Gust  case:  M=0.30  at  sea  level 


Semi-span  =  31.14  m 


Figure  A4:  Rectangular  wing  geometry 
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The  lift  distribution  was  computed  writing  the  previous  equations  in  Matlab  and  solving  the 
system.  The  result  is  obtained  by  dividing  the  wing  span  in  20  sections.  The  calculated  lift 
distribution  is  presented  in  the  Figure  A8. 


Figure  A5:  Lift  Distribution  of  the  rectangular  wing  model  (M=0.3,  AoA=0°) 

The  lift  distribution  obtained  with  this  theory  has  an  elliptical  shape.  This  lift  distribution  is 
compared  with  the  other  methods  further  in  the  chapter  where  the  effects  of  the  different 
assumptions  are  discussed. 


A.2.2  The  Weissinger  Theory 

The  Weissinger  Theory  [23]  is  based  on  the  Lifting  Line  Theory  but  accommodates  it  to  swept- 
back  wings.  Coefficients  taking  into  account  the  swept  back  angle  are  added  to  the  circulation 
equations  presented  in  the  Lanchester-Prandtl  Theory.  The  total  circulation  is  the  circulation  from 
the  Lifting  Line  Theory  plus  a  correction  term  which  takes  into  account  the  sweep  of  the  wing. 
The  ‘Influence  function’  L  is  used  for  the  correction  term: 


UKy0-y)] 


Vl  +  A2(y0-y)2  -  1 

K%-y) 


(A.ll) 


Where  : 


A(y)  =  is  the  local  aspect  ratio 

C(y) 


y0  =  57-  dimensionless  coordinate  in  the  absolute  plane 

/  2 
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So  the  equation  (A.3)  becomes: 


f  dy  [2  +  A(y0  -  y)L[A(y0  -  y)]\dy 

j-iyo  —  y 


However,  this  theory  is  more  complex  to  apply  than  the  Lifting  Line  Theory  and  would  be  time- 
consuming.  Thus,  the  author  preferred  to  use  directly  computational  methods  to  determine  and 
compare  the  lift  distribution  to  the  basic  theory  of  Prandtl. 

A.2.3  Computational  Analyses:  XFLR5 

The  aerodynamic  studies  were  performed  using  two  different  softwares,  FLUENT  CFD  and 
XFLR5.  In  parallel,  the  author  computed  the  same  cases  in  XFLR5. 

The  first  step  in  the  XFLR5  software  is  to  set  up  the  aerofoil  of  the  wing.  This  analysis  in  2D  is 
performed  using  the  software  X-foil,  included  in  the  XFLR5  package. 

The  aerofoil  NACA  4415  was  studied  for  different  conditions,  i.e.  various  Reynolds  numbers  and 
angles  of  attack.  The  Reynolds  number  depends  on  the  study  case  and  the  wing  chord. 


VL 


(A.13) 


Re  =  — 


v 


Where  L  is  a  characteristic  linear  dimension,  here  the  local  chord. 

The  range  of  the  Reynolds  number  for  the  geometry  is  very  large.  The  minimum  Reynolds 
number  is  located  from  the  kink  at  10.5  m  spanwise  to  the  wing  tip,  where  the  chord  is  the 


smallest  and  reaches  the  value  of  7.106  for  the  cruise  case.  The  maximum  value  of  the  Reynolds 


number  is  1.108  at  the  root  chord  for  the  gust  case.  For  all  cases,  the  Reynolds  number  is  high 
along  the  geometry.  The  viscosity  is  a  characteristic  of  the  fluid  to  consider  for  this  study  and  in 
particular  for  high  angles  of  attack  where  flow  separation  occurs.  The  X-foil  analysis  of  the  lift 
coefficient  for  the  NACA  4415  has  been  conducted  at  different  Reynolds  numbers.  The  Reynolds 
number  used  by  X-foil  takes  the  characteristic  dimension,  the  chord,  as  one.  Thus,  the  range  of 
Reynolds  numbers  to  be  calculated  by  X-foil  at  the  different  cases  is  between  1.6.106  and  7.106 


[24], 
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Figure  A6:  Lift  coefficient  vs.  a  (at  different  Reynolds  numbers,  by  unit  chord) 

The  software  XFLR5  computes  aerodynamic  data  in  three  dimensions.  Once  the  aerofoil  was 
defined  in  the  software,  the  wing  geometry  was  created  following  the  coordinates  given  in  the 
previous  paragraphs  for  the  new  geometry.  Several  analyses  were  conducted  to  determine  the 
spanwise  lift  distribution  and  pitching  moment.  The  local  lift  coefficient  CL  and  the  local  pitching 
moment  coefficient  Cm  were  obtained  using  both  the  Vortex  Lattice  Method  (VLM)  and  the  3D 
Panels  Method.  The  Vortex  Lattice  Method  and  3D  Panels  Method  are  two  methods  available  in 
the  software  to  compute  the  aerodynamic  loads. 


Figure  A7:  Wing  geometry  in  XFLR  5 


35 


Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


2013 


UNIVERSITY 


The  VLM  simplifies  the  wing  geometry  in  an  equivalent  2D  geometry  divided  into  panels.  A 


discrete  horseshoe  vortex  is  applied  at  the  control  point  of  each  panel.  The  theory  of  the  VLM  is 


based  on  the  Laplace’s  equations  and  is  an  extension  of  the  Lifting  Line  Theory.  For  each  panel, 
the  velocities  and  singularities  induced  by  the  vortex  are  computed.  Thus,  the  pressure  on  the 
surfaces,  the  lift  and  the  drag  are  calculated. 

However  this  method  does  not  take  into  account  the  viscosity  and  compressibility  of  the 
airflows.  These  assumptions  reduce  the  range  of  the  applications,  since  only  subsonic  flow  can  be 
modelled  (Mach  numbercl).  Although  the  VLM  relies  on  the  theory  of  ideal  flow  and  thus  the 
Laplace’s  equations,  the  compressibility  of  the  flows  can  be  corrected  for  high  subsonic  speeds  by 
using  the  Prandtl-Glauert  transformation: 


The  formula  can  be  used  for  the  lift,  drag  and  pitching  moment  coefficients  as  they  are  linear 
integrals  of  the  pressure  coefficient. 

On  the  other  hand,  since  the  viscosity  of  the  fluid  is  not  considered  in  the  calculations,  the  skin- 
friction  drag  is  not  added  to  the  total  drag.  The  influence  of  the  thickness  is  not  taken  into  account 
in  the  calculations  using  VLM  as  well.  However,  the  3D  Panels  Method  models  the  wing  in  three 
dimensions,  discretizing  the  span  into  panels  following  the  aerofoil  curve.  The  curves  are 
idealized  by  straight  lines  at  each  aerofoil  section.  Therefore,  both  upper  and  lower  surface 
characteristics  are  considered  separately  in  the  calculations.  Although  the  3D  Panels  Method  has 
the  same  theoretical  restrictions  as  the  Vortex  Lattice  Method,  this  method  has  the  advantage  of 
considering  the  thickness  of  the  aerofoil. 

For  low  Mach  numbers  (M<0.3)  and  high  Reynolds  numbers  (Re>105),  these  assumptions  can 
be  done  to  obtain  initial  results  using  these  relative  easy  and  rapid  methods  [25].  The  lift 
distribution,  pitching  moment  distribution  and  location  of  the  centre  of  pressure  along  the  span 
were  extracted  from  the  two  methods.  The  results  are  compared  to  the  other  methods,  and 
particularly  to  the  CFD  results  in  the  next  paragraph.  However,  additional  calculations  were 
needed  after  the  acquisition  of  the  data  from  the  software  in  order  to  obtain  the  desired 
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information.  For  example,  the  location  of  the  aerodynamic  centre  along  the  span,  which  is 
important  for  the  load  calculation  (cf.  Appendix  B),  was  calculated  from  the  results  of  the  centre 
of  pressure. 


Section 

Aerofoil  Approximation 


Figure  A8:  3D  Panels  approximation  to  an  aerofoil 


From  the  two  methods  used  in  XFLR5,  the  centre  of  pressure  was  located  along  the  span.  The 
centre  of  pressure  is  the  point  where  the  total  lift  of  the  section  applies  without  creating  any 
moment. 

Although  the  aerodynamic  centre  is  often  assumed  at  one  quarter  of  the  chord,  the  aerodynamic 
centre  location  for  this  geometry  was  computed  from  the  centre  of  pressure  location  for  different 
angles  of  attack  in  order  to  verify  this  common  assumption. 

The  aerodynamic  centre  is  the  point  where  the  pitching  moment  coefficient  of  the  section  does  not 
vary  with  the  lift  coefficient: 


d4r  = 0  <A15> 

dCL 

Since  the  lift  coefficient  at  a  given  station  depends  only  on  the  angle  of  attack  at  given  flow 
conditions,  and: 


L'  =^pcCLV2 
M'  =\pc2CmV2 


(A.16) 

(A.17) 
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Figure  A9:  Moment  at  the  aerodynamic  centre 


As  presented  in  the  Figure  A 12,  the  aerodynamic  centre  location  can  be  calculated  as  following: 


L'l(.XcPl  —  XAC )  —  L' 2  (,XCP2  —  XAC ) 


(A.18) 


So, 


XAC  — 


^2XCP2  ^jAcPl 

l' 2  - 


(A.19) 


Once  the  results  were  obtained  with  XFLR5,  the  lift  distribution,  pitching  moment  distribution, 
aerodynamic  centre  location  were  then  compared  to  the  CFD  analyses  results. 

A.3  Comparison  to  the  CFD  Analysis 

The  VLM,  Lifting  line  and  3D  Panels  method  results  can  be  compared  to  available  in  house 
Computational  Fluid  Dynamics  (CFD)  results  which  used  the  same  geometry.  The  commercial 
software  Fluent  was  used.  The  calculations  with  CFD  take  into  account  the  viscosity  and 
compressibility.  The  mesh  can  be  denser  at  critical  points  in  the  flow  field  to  capture  the  flow 
physics  better.  This  is  especially  important  for  regions  with  high  gradients. 

The  results  of  the  different  analysis  methods  described  previously  were  compared  for  the  two 
study  cases.  The  initial  data  gave  an  angle  of  attack  of  0°  for  both  cases. 

A.3.1  Gust  Case:  Sea  Level,  Mach  0.3 

The  Figure  A13  presents  the  results  obtained  for  the  spanwise  lift  distribution  using  different 
methods  for  the  gust  case  with  an  angle  of  attack  of  0°.  The  VLM  and  3D  Panels  method  give 
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values  very  close  to  the  CFD  calculations  with  a  maximum  difference  of  10%.  The  shape  of  the 
distributions  is  quite  the  same  except  at  the  tip  of  the  wing  where  we  can  observe  a  drop  of  the  lift 
at  31  m  of  the  span  in  the  CFD  analysis  whereas  the  VLM,  3D  Panels  lift  distributions  keep  a 
smooth  decreasing  curve. 


Figure  A10:  Spanwise  lift  distribution  for  gust  case  (AoA=0°) 


The  Lifting  Line  Theory  gives  values  close  to  the  other  methods  for  the  straight  swept  part 
from  10.5  m  to  the  tip  even  if  this  method  tends  to  simplify  the  lift  distribution  by  doing  many 
assumptions.  However,  the  difference  in  the  accuracy  of  the  method  compared  to  the  others  is 
highlighted  by  under  evaluating  the  lift,  reaching  a  difference  of  20%  at  the  kink  (10.5  m 
spanwise).  The  results  from  the  root  to  the  kink  are  not  relevant  with  this  method  and  cannot 
compared  to  the  other  method’s  results  as  the  idealized  geometry  presented  in  the  paragraph  A.2.1 
represents  only  the  part  of  the  geometry  from  the  kink  to  the  tip  of  the  wing. 

As  the  total  lift  created  by  the  wing  was  higher  than  the  weight  that  balanced  the  aircraft,  the  angle 
of  attack  of  the  aerofoil  was  reduced  to  -1°.  With  this  new  configuration,  the  equilibrium  is 
maintained. 
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Figure  A1 1 :  Spanwise  lift  distribution  for  gust  case  (AoA=-1  °) 

The  new  lift  distributions  have  the  same  shape  as  the  previous  one  for  zero  angle  of  attack.  The 
total  lift  has  been  reduced  from  380  kN  to  300  kN. 


Figure  A12:  Location  of  the  aerodynamic  centre  for  gust  case 

It  can  be  highlighted  in  the  Figure  15  that  all  the  methods  mostly  calculate  the  aerodynamic  centre 
located  between  24%  and  28%  of  the  chord  along  the  span  except  at  the  tip  where  the  tip  vortex 
affects  the  lift  distribution.  The  3D  Panels  method  diverges  from  the  other  results  at  the  root  but 
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follows  the  same  trend  as  the  CFD  at  the  wing  tip.  From  the  CFD  results,  it  can  be  showed  that  the 
aerodynamic  centre  moves  backward  at  the  tip,  located  at  39%  of  the  chord.  This  phenomenon  has 
already  noticed  in  previous  experiments  [20],  Nevertheless,  the  common  assumption  which  gives 
the  location  of  the  aerodynamic  centre  at  the  quarter  chord  can  be  validated  for  most  of  the  span. 
Therefore,  the  exact  values  of  the  aerodynamic  centre  were  taken  from  the  CFD  for  the  following 
studies,  as  the  CFD  is  considered  as  the  most  accurate  method. 


Spanwise  Pitching  Moment  Distribution  AoA=-l° 


Semi-Span  (m) 


3D  panels 
CFD 


Figure  A13:  Spanwise  pitching  moment  distribution  for  gust  case  (AoA=-1  °) 

The  three  methods  CFD,  VLM  and  3D  panels  give  the  same  results  of  pitching  moment  from  10.5 
m  to  the  tip  of  the  wing  which  corresponds  to  the  straight  swept  part  of  the  wing.  The  pitching 
moment  decreases  rapidly  to  the  root  of  the  wing.  However,  the  VLM  results  diverge  from  the 
other  results  for  the  triangular  part  of  the  wing. 


A.3.2  Cruise  Case:  60000  ft,  Mach  0.65 

As  done  for  the  gust  case,  the  lift  distribution  was  computed  for  the  cruise  case  (altitude  60000  ft 
and  M=0.65).  The  CFD,  VLM,  3D  Panels  and  Lifting  Line  methods  results  are  compared  in  the 
Figure  17  for  the  angle  of  attack  0°  as  given  in  the  specifications. 
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Figure  A14:  Spanwise  lift  distribution  for  cruise  case  (AoA=0°) 

The  results  from  the  different  methods  are  close  and  follow  the  same  trends.  As  explained  for 
the  gust  case,  only  the  results  from  10.5  m  spanwise  are  relevant  for  the  lifting  line  theory.  The 
discrepancy  between  the  wing  tip  lift  calculated  using  XFLR5  and  CFD  increases  compared  to  the 
gust  case.  The  discrepancy  reaches  14%  between  the  3D  Panels  results  and  the  CFD  results  at 
28.5  m  spanwise. 

However,  the  analyses  revealed  that  the  total  lift  created  by  the  wing  in  these  conditions  do  not 
reach  a  sufficient  value  to  counterbalance  the  weight  of  the  aircraft.  In  these  conditions,  the  lift 
decreases  by  more  than  half  compared  to  the  gust  case.  The  total  lift  is  131  kN  which  is  not 
enough  for  the  maximum  weight  of  271  kN,  even  if  it  is  assumed  that  20  %  of  the  fuel  has  been 
burned  to  reach  this  altitude. 

In  order  to  achieve  the  lift  needed  for  the  cruise,  analyses  have  been  computed  for  the  same 
freestream  conditions  but  increasing  the  angle  of  attack.  As  presented  in  the  Figure  A9,  the  lift 
coefficient  CL  increases  when  the  angle  of  attack  increases  between  zero  and  ten  degrees. 

The  analyses  have  shown  that  an  angle  of  attack  of  5°  is  required  to  have  enough  lift  for  the  case. 
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Spanwise  Lift  Distribution  AoA=5° 
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Figure  A15:  Spanwise  lift  distribution  for  cruise  case  (AoA=5°) 


It  can  be  noticed  in  the  Figure  30  that  the  lift  distribution  calculated  from  XFLR5’s  methods  and 
the  Lifting  Line  Theory  diverge  from  the  CFD  results.  This  is  especially  visible  for  the  straight 
swept  part  from  10.5  m  to  the  tip.  The  lower  lift  from  the  CFD  can  be  attributed  to  a  local  flow 


separation  at  the  kink  in  the  geometry. 


Figure  A16:  Location  of  the  aerodynamic  centre  for  cruise  case 


Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


43 


2013 


EOARD/AFRLl 


Cranfield 

I  UNIVERSITY 


As  for  the  gust  case,  the  aerodynamic  centre  is  located  very  close  to  the  quarter  of  the  chord 
(cf.  Figure  A 19).  The  location  of  the  aerodynamic  centre  moves  aft  at  the  tip  of  the  wing.  The 
CFD  gives  the  aerodynamic  centre  at  34%  of  the  chord,  which  is  due  to  the  tip  vortex. 

The  pitching  moment  shape  does  not  change  between  the  methods.  However,  as  for  the  lift 
distribution,  the  local  pitching  moment  value  decreased  by  half  compared  to  the  gust  case  (cf. 
Figure  A20).  This  is  explained  by  the  decrease  of  the  dynamic  pressure.  The  density  of  the  air 
drops  from  1.225  kg/m3  at  sea  level  to  0.115  kg/m3  at  60000  ft  when  velocity  is  just  doubled. 


Figure  A17:  Spanwise  pitching  moment  distribution  for  cruise  case  (AoA=0°) 

This  chapter  has  compared  the  results  for  the  lift,  pitching  moment  distribution  and  the 
aerodynamic  centre.  Despite  of  the  assumptions  done  by  the  VLM  and  3D  panels  method,  the 
results  have  shown  convergence  with  the  CFD  results. 

The  CFD  results  have  been  used  for  the  following  load  calculations.  It  would  be  beneficial  study 
to  perform  wind  tunnel  tests  to  validate  the  computational  results. 
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Appendix  B.  Mass  Estimation  and  Load  Calculation 

B.l  Mass  Distribution 

The  mass  estimation  of  the  wing  consisted  in  three  distinct  parts:  the  structural  mass,  the 
system  mass  and  the  fuel  mass.  From  the  specifications,  no  detailed  data  was  given  on  the 
proportion  of  the  structural  mass  over  the  total  mass  or  on  the  location  of  the  systems.  Only  the 
fuel  mass  was  known.  The  author  decided  to  make  assumptions  on  the  different  masses  and 
locations  of  the  systems  as  well  as  on  the  structural  mass.  The  procedure  adopted  to  evaluate  these 
parameters  is  explained  in  the  next  chapters. 

B.1.1  Structural  Mass  Distribution 

As  the  system  mass  distribution  was  not  given,  the  easiest  way  to  evaluate  the  structural  and 
system  masses  was  to  evaluate  the  structural  mass  first. 

The  wing  model  was  approximated  as  a  rectangular  beam  representing  both  front  and  rear  spars, 
on  which  the  lift  is  applied.  The  beam  was  sized  to  support  the  bending  due  to  the  lift  at  ultimate 


loads. 


The  height  of  the  beam  is  the  average  value  of  the  front  and  rear  spar  heights  located  respectively 
at  15%  and  75%  of  the  chord  at  the  tip  of  the  wing  and  kept  parallel  to  the  leading  edge.  The 
layout  of  the  wing  structure  is  described  in  further  details  in  the  next  chapter. 


h 


< — > 


b 


Figure  18:  The  beam  under  the  bending  moment 
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The  thickness  of  the  beam  has  to  be  large  enough  not  to  fail  under  loading.  The  ultimate  stress  that 
the  beam  can  support  is  described  by  the  following  criterion: 


_  Mh 
°uit  —  2i 


(B.2) 


Where: 


I  =  —  is  the  second  moment  of  inertia  of  the  rectangular  beam 
M  is  the  ultimate  bending  moment 

Gun  =  500  MPa  is  the  ultimate  strength  of  composite  materials,  taken  as  an  average  value 
Section  by  section,  the  thickness  was  computed  and  the  area  of  the  beam  was  estimated.  The  mass 
of  the  beam  was  calculated  using  a  density  of  composite  materials  of  1600  kg/m3. 


Structural  Mass  Distribution  Estimation 


Figure  19:  Estimation  of  the  structural  mass  distribution 

The  drop  of  the  mass  per  unit  of  span  at  10.5  m  spanwise,  i.e.  at  the  kink,  is  due  to  the  significant 
increase  of  the  height  of  the  wing  box  which,  consequently,  has  a  higher  moment  of  inertia  I. 

By  integrating  this  mass  distribution  along  the  wing  span,  a  total  mass  of  5700  kg  was  estimated 
for  the  wing  structure. 

B.1.2  System  Mass  Distribution 

The  system  mass  was  deduced  from  the  structure  estimation,  by  subtracting  it  from  the  total 
empty  mass. 
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Comparing  to  other  aircraft,  the  engine  mass  was  assumed  to  be  2  tons  per  engine,  with  one 
engine  per  side.  The  engine  is  located  at  5  m  outboard  from  the  centreline. 

The  landing  gears  weigh  2.2  tons  in  total  (both  sides).  The  main  landing  gear  is  considered  to 
represent  90%  of  the  mass  and  located  at  9m  from  the  nose.  The  nose  landing  gear  is  at  1.5  m  aft 
from  the  datum  on  the  centreline.  The  location  of  the  engines  and  landing  gears  are  presented  in 
the  Figure  24. 


y 


Figure  20:  Layout  of  the  main  components  and  fuel  tanks 


The  remaining  systems  mass  consists  in  the  actuation  system,  fuel  and  oil  pumps,  auxiliary  power 
unit,  avionics  systems... 

20%  of  this  mass  was  distributed  from  the  engine  bays  to  the  wing  tips  for  the  actuation  and  fuel 
system.  The  other  80%  was  distributed  inboard,  between  the  engine  bays.  The  mass  was 
distributed  proportionally  to  the  volume  of  the  cross-sections  of  the  wing. 
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Systems  Mass  Distribution 


Figure  21 :  Estimation  of  the  system  mass  distribution 
B.1.3  Fuel  Mass  Distribution 

Several  fuel  tanks  were  assumed  along  the  span  to  fulfill  the  total  fuel  mass  of  16  tons  (for  half 
of  the  aircraft).  The  inboard  fuel  tank  1  (in  red  in  the  Figure  24)  is  delimited  by  the  avionic  bay 
spar.  The  fuel  tanks  2  and  3  are  restricted  between  the  front  and  rear  spars.  This  configuration  lets 
enough  space  for  the  systems  inboard  and  permits  to  have  a  maximum  of  load  of  fuel  forward. 
This  choice  has  been  done  also  to  have  the  maximum  of  fuel  in  the  front  part  of  the  inboard  wing 
to  move  the  centre  of  gravity  forward  for  stability  reasons.  The  outboard  tank  3  is  stopped  3  m 
before  the  tip  of  the  wing  to  allow  space  for  the  gust  alleviation  device. 

The  inboard  fuel  tank  volume  is  large  enough  to  contain  almost  all  the  fuel.  Having  a 
maximum  of  fuel  inboard  is  not  the  best  solution  to  reduce  the  loads  and  bending  of  the  wing.  The 
author  preferred  to  distribute  the  fuel  more  uniformly  along  the  span.  As  a  consequence,  the 
inboard  fuel  tanks  are  filled  at  only  25%  of  their  capacity.  The  outboard  tanks  are  full  of  fuel,  with 
20%  of  the  volume  measured  in  CATIA  reserved  to  take  into  account  the  structure  and  systems. 
The  final  fuel  and  total  mass  distributions  are  presented  in  the  Figure  26. 
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Wing  Mass  Distribution 


Figure  22:  Estimation  of  the  system  mass  distribution 

B.2  Load  Calculation 

B.2.1  Aerodynamic  Data 

The  data  obtained  from  the  CFD  was  used  to  compute  the  loads  on  the  wing.  When  CFD 
results  were  not  available,  3D  panels  results  were  used  instead.  This  is  the  case  for  the  gust  case  at 
empty  weight  where  an  angle  of  attack  of  -3°  is  necessary  to  balance  the  weight.  The  second  sea 
level  case  needs  an  angle  of  attack  of  -1°  from  the  CFD  results. 

For  the  cruise  cases,  the  angles  of  attack  of  5°  and  0°  has  been  respectively  chosen  for  full-fuel 
weight  and  the  empty  weight  cases  to  produce  enough  lift  to  compensate  the  weight. 

In  order  to  keep  the  equilibrium  of  the  aircraft,  the  lift  distributions  were  scaled  down  to  obtain 
the  exact  amount  of  lift  needed  to  compensate  the  weight. 

In  real  conditions,  the  angle  of  attack  would  be  adapted  to  the  conditions  to  compensate  exactly 
the  weight  and  obtain  steady  level  flight.  Furthermore,  the  lift  might  be  reduced  by  the  deflection 
of  the  elevator  and  other  surfaces.  The  aim  of  the  aerodynamic  analyses  was  to  provide  the  lift 
distribution  shape  and  confirm  that  enough  lift  can  be  produced  by  the  wing  for  the  different 
cases. 
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- Gust  case  1 

(AoA=-l°) 

- Gust  Case  2 
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—  •  -  Cruise  case  1 
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Figure  23:  Lift  distributions  used  in  the  loading  actions 


B.2.2  Load  Factors 

The  cases  studied  in  this  analysis  occur  at  different  flight  conditions.  The  maximum  load  factor 
the  structure  can  encounter  in  each  situation  needs  to  be  calculated  to  determine  the  worst  loads. 
With  regards  to  the  span  of  the  aircraft  and  the  take-off  mass,  the  certifications  used  for  the 
analyses  are  the  EASA  CS-25,  the  certifications  for  large  aeroplanes. 

First  of  all,  the  flight  envelope  has  to  be  defined.  The  maximum  load  factor  n  of  the  flight 
envelope  that  the  structure  has  to  support  is  calculated  from  the  equations  below: 

n>  2.1 +  ( — — — )  (MTOW  in  lb)  (B.1) 

VMrow+ioooo/ '  '  '  ' 

And  2.5  <  n  <  3.8 

This  gives  a  maximum  load  factor  of  2.5  for  the  flight  envelope. 

For  each  case,  the  gust  load  factor  is  determined  for  initial  structural  design.  The  gust  load  factor 
was  calculated  from  the  alleviated  sharp  edge  analysis  presented  in  the  CS-23. 

From  the  article  CS-23. 333,  the  gust  velocities  Ude  are  for  each  case: 

Gust  case  (at  sea  level):  Ude  =  50  ft/s 
Cruise  case  (at  60000  ft):  Ude  =  25  ft/s 
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The  lift  curve  slope  ax  was  deduced  from  the  aerofoil  aerodynamic  study  done  with  X-Foil  and 
confirmed  with  the  NACA  Report  N°832  [26]. 

The  results  are  presented  in  the  next  table: 


Table  1 :  Gust  load  factor  for  the  study  cases 


Gust  case  1 

Gust  case  2 

Cruise  case  1 

Cruise  case  2 

Gust  load  factor  2.95 

3.81 

2.13 

3.42 

In  order  to  calculate  the  maximum  loads  that  the  structure  has  to  support,  the  final  load  factor  is 
chosen  between  the  gust  load  factor  and  the  flight  envelope  load  factor  of  2.5  whichever  is  the 
greatest. 

B.2.3  Shear  Force,  Bending  Moment,  Torque  Diagrams 

Once  the  lift  and  mass  distributions  were  obtained,  the  shear  force  and  bending  moment  were 

calculated  along  the  span. 

First  of  all,  the  convention  used  in  the  diagrams  is  illustrated  in  Figure  27. 


Figure  24:  Positive  convention  for  the  shear  force,  bending  moment  and  torque  diagrams 

The  shear  force  was  calculated  section  by  section,  starting  by  the  tip  of  the  wing  where  the 
loads  are  equal  to  zero.  Shear  force  distributions  were  first  computed  for  a  load  factor  of  lg. 
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Figure  25:  Shear  force  diagram  for  the  gust  case  2  (1  g) 

Then,  the  results  were  multiplied  by  the  load  factors  calculated  previously  for  the  different 
flight  conditions.  In  the  Figure  29,  the  shear  force  diagrams  for  the  different  cases  are  plotted  at 
limit  load,  considering  the  individual  load  factors.  The  shear  force  starts  increasing  from  the  tip 
where  the  lift  is  locally  higher  that  the  weight.  It  can  be  highlighted  that  the  shear  force  drops 
suddenly  at  5m  for  all  the  cases.  This  is  due  to  the  weight  of  the  engine  which  is  predominant  over 
the  lift.  Then,  the  fuel  weight  and  the  landing  gear  weight  make  the  shear  force  decrease  to  zero  at 
the  centreline,  which  is  expected  as  the  total  lift  is  equal  to  the  total  weight.  The  gust  case  2,  the 
gust  case  at  sea  level  with  no  fuel,  gives  the  highest  shear  force  for  a  large  part  of  the  span. 
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Shear  Force 


Figure  26:  Shear  force  diagram  for  the  different  cases  (limit  loads) 

The  bending  moment  is  the  integration  section  by  section  of  the  shear  force  from  the  tip  to  the 
root  along  the  elastic  axis.  The  torque  is  also  calculated  on  the  elastic  axis,  assumed  to  be  located 
at  40%  of  the  chord.  The  moment  at  the  shear  centre  is  obtained  from  the  aerodynamic  load  and 
moment  at  the  aerodynamic  centre  and  the  weight  at  the  centre  of  gravity  of  each  section.  As  the 
wing  is  swept,  corrections  on  the  bending  and  torque  have  to  be  done  to  obtain  it  in  the  local  axis, 
aligned  with  the  elastic  axis  [27]. 

The  diagrams  for  different  cases  are  represented  in  the  Figure  30  and  31.  The  bending  increases 
gradually  to  reach  its  maximum  value  at  the  centreline.  The  gust  case  2  is  again  the  worst  case 
with  a  maximum  bending  of  4400  kN.m  at  the  centreline.  The  torque  is  positive  (nose  up)  from 
the  kink  at  10.5  m  to  the  tip  because  of  the  correction  done  to  have  the  values  in  the  local  axis. 
Indeed,  corrected  torque  is  calculated  with  an  equation  taking  into  account  both  non-corrected 
torque  and  bending  moment  values  and  the  sweep  angle  of  the  elastic  axis.  The  bending  moment 
is  predominant  in  the  equation  for  this  portion  so  that  the  final  torque  becomes  positive.  Then, 
from  the  kink  to  the  centreline,  the  pitching  moment  increases  significantly  (cf.  Appendix  A)  and 
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creates  a  nose-down  torque.  The  fuel  tanks  located  upstream  in  the  inboard  part  of  the  wing 
amplify  this  effect  as  well.  The  effect  of  the  aft  location  on  the  chord  of  weight  of  the  engine  and 
main  landing  gear  is  observed  at  4  m  and  5  m  of  the  span  by  reducing  the  nose-down  torque. 


Figure  27:  Bending  moment  diagram  for  the  different  cases  (limit  loads) 


Figure  28:  Torque  diagram  for  the  different  cases  (limit  loads) 

The  envelopes  of  the  different  diagrams  were  used  for  the  initial  sizing  of  the  structure. 
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Appendix  C.  Initial  Structural  Layout 

C.l  General  Layout 

As  briefly  presented  in  the  mass  distribution  chapter,  the  preliminary  design  of  the  flying  wing 
layout  has  been  carried  out.  The  structure  aims  to  support  the  loads  by  creating  major  load  paths 
through  the  skins,  spars,  ribs  and  frames.  The  total  structural  weight  depends  on  the  configuration 
and  the  arrangement  of  the  different  members. 


Figure  29:  Main  components  of  the  structure 

The  wing  can  be  divided  into  two  distinct  parts  considering  its  geometry:  the  inboard  wing  and  the 
outboard  wing.  The  inboard  wing  is  constituted  of  the  large  tapered  central  part  going  to  the  kink 
at  around  10.5  m  laterally.  The  outboard  wing  is  the  straight  swept  back  part  from  the  kink  to  the 
tip. 
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C.2  Outboard  Wing 

The  outboard  wing  has  a  conventional  two-spar  configuration.  The  front  spar  and  the  rear  spar  are 
respectively  located  at  15%  and  75%  of  the  chord  at  the  wing  tip  and  are  kept  parallel  to  the 


leading  edge  /  trailing  edge  along  the  span.  The  upper  and  lower  skin  are  reinforced  by  I-section 


stringers  to  avoid  the  buckling  as  high  bending  moments  are  applied  to  the  covers  creating  high 
compressive  loads.  For  the  upper  skin,  subjected  to  the  highest  compressive  loads,  the  stringer 
pitch  has  been  set  to  200  mm,  following  the  typical  values  in  aircraft  design  [13].  Contrary  to  the 
upper  skin,  the  compressive  loads  are  lower  for  the  bottom  skin  as  the  design  case  is  the  bending 
due  to  the  weight  of  the  wing  (lg  case).  In  the  objective  of  weight  reduction,  the  stringer  pitch  has 
been  increased  to  400  mm  for  the  lower  skin  from  the  section  8  to  the  tip  (cf.  Figure  34).  A 
balance  between  the  number  of  stringers  and  the  skin  thickness  has  been  defined  to  have 
proportional  dimensions  between  the  skin  and  stringers  designs.  The  skin  also  needs  a  minimum 
thickness  regarding  lightning  protection.  The  ribs  are  perpendicular  to  the  rear  spar  to  give  the 
best  arrangement  in  terms  of  load  transmission  from  the  trailing  edge  devices.  The  rib  pitch  is  set 
at  750  mm  at  the  kink  and  increased  progressively  to  1  m  at  the  tip  of  the  wing,  where  the  loads 
are  lower. 

The  tip  of  the  wing  has  been  considered  separately  in  the  design  phase  to  take  into  account  the 
future  presence  of  the  passive  gust  alleviation  device.  The  last  three  meters  are  reserved  for  the 
device  design.  The  front  and  rear  spars  are  kept  at  the  same  emplacement  but  the  rib  direction  is 
changed  for  the  streamline  direction.  This  choice  has  been  motivated  by  the  absence  of  trailing 
edge  devices  and  the  fact  that  the  device  rotates  around  the  shaft  in  the  y-axis.  As  the  bending  and 
torque  at  the  tip  are  low,  the  stringer  pitch  has  been  increased  to  400  mm  for  both  upper  and  lower 
skin. 
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Figure  30:  Upper  skin  layout 


C.3  Inboard  Wing 

The  inboard  wing  is  composed  of  numerous  elements.  With  its  triangular  shape,  a  conventional 
wing  design  cannot  be  applied  in  this  case.  The  loads  are  coming  for  the  outboard  part  of  the  wing 
through  the  front  and  rear  spars  and  the  stiffened  skin  panels. 

In  order  to  keep  continuity  of  the  load  paths,  the  stringers  follow  the  same  direction  and  shape  as 
the  inboard  part.  Since  the  height  of  the  wing  box  increases  significantly  from  the  kink  to  the 
centreline  (0.5  m  at  the  kink  to  1.5  m  at  the  centreline)  and  the  applied  compressive  loads  are 
inversely  proportional  to  the  wing  box  height,  the  skin  panels  do  not  need  to  be  as  stiff  as  the 
outboard  skins.  Thus,  the  stringer  pitch  has  been  increased  to  400  mm  for  both  upper  and  lower 
covers,  reducing  the  weight  of  the  covers.  The  front  and  rear  spars  of  the  outboard  wing  come 
through  the  fore  part  of  the  wing.  Additional  spars  called  ‘middle  spars’,  are  located  in  the  centre 
of  the  wing.  They  separate  the  inboard  wing,  delimit  the  fuel  tanks  and  help  supporting  the 
spanwise  shear  force,  the  engine  and  landing  gear  masses  (cf.  Figure  24).  A  ‘second  rear  spar’ 
comes  from  the  kink  and  follows  the  trailing  edge  to  give  a  support  for  the  control  devices  and 
flaps.  Frames  are  designed  to  support  the  chordwise  bending.  The  heavy  frames  are  constituted  of 
a  web  and  I-section  caps  attached  to  the  upper  and  lower  skins.  They  support  the  chordwise 
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bending  as  well  as  the  chordwise  shear.  They  have  the  same  function  as  the  ribs  and  help 
supporting  the  heavy  mass  of  the  systems  such  as  the  engines.  The  light  frames,  set  as  chordwise 
stiffeners,  are  placed  between  the  heavy  frames  to  protect  the  skin  against  buckling.  The  light 
frames  are  I-section  beams  uniformly  spaced  of  approximately  750  mm  between  the  heavy 
frames.  They  are  boundaries  for  the  spanwise  stringers. 
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Appendix  D.  Initial  Sizing 

D.l  Material  Selection 

The  selection  of  the  material  is  an  important  decision  to  make  as  it  has  a  direct  impact  on  the 
structural  design  and  stress  analysis  and  thus  on  the  weight  of  the  structure.  Weight  saving  is  the 
first  concern  of  all  aircraft  designers  to  reduce  the  fuel  consumption,  increase  the  maximum  range 
and  the  payload.  Aluminium  alloys  have  been  the  most  chosen  for  several  decades;  their 
properties  have  been  improved  over  time.  However,  the  composite  materials  take  an  increasingly 
important  role  in  the  structure  materials. 

The  composite  material’s  reduced  density  is  the  one  of  the  advantages  of  these  materials 
compared  to  the  aluminium  alloys  but  their  anisotropic  properties  make  the  design  complex.  For  a 
flying  wing,  the  high  tensile  performance  of  the  CFRP  laminates  is  ideal  for  the  skin  design  which 
is  subject  to  high  bending.  The  laminate  layup  has  to  be  tailored  to  the  main  load  paths  of  each 
component  to  maximize  the  benefits  of  the  use  of  the  composite  materials.  The  bending-torsion 
coupling  properties  of  the  laminate  can  be  adapted  to  improve  the  aeroelastic  behaviour,  the 
structural  modes  of  the  wing.  The  composites  present  also  a  longer  fatigue  life  which  allows  either 
to  delay  the  regular  inspections  or  to  reduce  their  number. 

The  common  8552  epoxy  matrix  IM7  UD  carbon  fibre,  an  intermediate  modulus  carbon  fibre 
epoxy  matrix  composite  has  been  chosen  for  the  whole  flying  wing.  The  properties  of  the  material 
are  presented  in  the  Table  4. 


Table  2:  8552/IM7  Material  properties  [28] 


0°  tensile  modulus 

GPa 

Pi 

164 

90°  tensile  modulus 

GPa 

e2 

12 

Shear  modulus 

GPa 

G12 

5.8 

0°  tensile  strength 

MPa 

X, 

2724 

90°  tensile  strength 

MPa 

Y, 

111 

Shear  strength 

MPa 

S 

120 

0°  compressive  strength 

MPa 

Xc 

1690 
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90°  compressive  strength 

MPa 

Yc 

250 

Poisson  ratio 

- 

Vl2 

0.3 

Density 

g/cm3 

P 

1.57 

Ply  thickness 

mm 

t 

0.131 

Fibre  volume 

% 

Vf 

57.7 

Once  the  material  was  chosen,  the  initial  sizing  of  the  different  wing  components  was  fulfilled. 
These  component  sizes  were  then  input  in  the  Finite  Element  program  to  analyze  the  structure. 

D.2  Member  Initial  Sizing 

D.2.1  Introduction 

The  member  initial  sizing  was  derived  from  the  shear  force,  bending  moment  and  torque 
diagrams.  The  worst  loads  from  the  four  cases  studied  in  the  Chapter  4  were  used  to  calculate  the 
local  stresses  in  the  structure  members.  The  geometric  data  needed  for  the  calculations,  such  as 
the  wing  box  height,  the  rib  pitch,  were  obtained  from  CATIA. 

The  laminate  engineering  properties  were  derived  from  CoALA  [29]  ,  an  in-home  software,  in 
which  the  laminate  plies  were  input.  The  stiffness  matrices  are  also  computed.  This  software 
calculates  the  failure  indices  and  strains  of  each  ply  of  the  laminate  for  a  given  load.  The  failure 
indices  were  checked  to  be  below  one  under  ultimate  loads  (limit  load  x  1.5).  The  strains  were 
kept  under  3500  ps  at  limit  loads  for  damage  tolerance. 

The  laminate  layups  chosen  are  all  balanced  and  symmetric.  In  order  to  keep  an  acceptable 
number  of  layup  combinations,  only  0°,  90°  and  +457-45°  angles  were  used  for  the  ply  direction. 
Typical  laminate  layups  are  defined  for  the  components.  They  are  repeated  to  obtain  the  desired 
thickness  of  the  components. 

The  laminate  design  followed  the  usual  guidelines: 

A  minimum  of  10%  of  plies  in  each  of  the  four  directions 
No  0°  and  90°  consecutive  plies 

A  maximum  of  four  consecutive  plies  of  the  same  direction 
+457-45°  directions  used  for  the  exterior  plies  for  damage  tolerance 
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In  order  to  facilitate  the  sizing,  design  sections  have  been  created  along  the  wing  span.  The 
following  initial  sizing  of  the  components  is  based  on  this  section  division. 


Figure  31 :  Sections  for  the  initial  sizing 

D.2.2  Skin  /  Stringers 

The  upper  and  lower  surfaces  are  subject  to  high  spanwise  bending  moments.  Consequently,  the 
upper  skin  is  reinforced  by  I-stringers  to  support  the  compressive  stresses. 

In  order  to  simply  the  initial  sizing  analysis,  a  single  laminate  layup  [+45/02/-45/90]s  was  used  for 
both  skin  and  stringers.  0°  plies  need  to  be  in  major  proportion  in  the  skin  layup  to  support  the 
axial  loads  but  also,  a  non-negligible  percentage  of  +45°/-45°  plies  is  needed  for  the  shear  induced 
by  the  torque  around  the  wing  box. 

The  stiffened  panels  were  designed  to  resist  buckling  and  keep  the  strains  under  the  limit  of  3500 
ps.  The  ultimate  bending  moment  obtained  from  the  loading  actions  was  used  for  the  upper 
surface  buckling  design  and  the  limit  loads  for  the  strain  design  for  both  upper  and  lower  surfaces. 
The  lower  surface  is  in  compression  when  the  aircraft  is  on  the  ground.  Thus,  the  bending  due  to 
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the  wing  weight  (lg  load)  was  computed  to  determinate  the  compression  loads  and  to  size  the 
lower  stiffened  surface. 

The  load  applied  to  the  skin/stringers  for  each  section  is  determined  by: 

M 

Nx=—  (D.l) 

x  wh 

Where  Nx  is  the  edge  load  applied  to  the  skin/stringers  panel,  M  is  the  bending  moment,  h  is  the 
height  of  the  wing  box  and  w  is  the  width  of  the  wing  box. 

This  analysis  was  conducted  by  V.  Fu  who  determined  the  best  design  for  the  skin  and  the 
stringers  at  each  section  using  an  in-house  code. 


Table  3:  Skin  panels  thicknesses 


Section 

Upper  skin 
thickness  (mm) 

Lower  skin 
thickness  (mm) 

Section 

Upper  skin 
thickness  (mm) 

Lower  skin 

thickness 

(mm) 

1 

4.45 

3.67 

7 

5.24 

3.41 

2 

5.24 

2.88 

8 

4.72 

3.93 

3 

6.03 

3.14 

9 

4.19 

3.41 

4 

7.60 

5.24 

10 

3.14 

2.62 

5 

6.29 

4.45 

11 

2.10 

2.92 

6 

5.76 

3.93 

D.2.3  Spars 

The  spar  web  is  designed  to  support  shear.  From  the  beginning,  the  layup  [+45/0/-45/90]S  was 
chosen  for  its  fifty  per  cent  of  +45/-45  plies  resisting  the  shear.  The  shear  is  calculated  from  the 
worst  values  of  shear  force  and  torque  in  the  wing  box.  The  spar  configuration  is  quite  complex  in 
the  inboard  part  of  the  wing  where  there  are  up  to  seven  spars  at  the  section.  The  initial  sizing  of 
the  spar  web  was  first  achieved  for  the  outboard  part  where  a  simple  front  and  rear  spar 
configuration  exists. 

The  following  equations  were  used  to  determine  the  shear  flow  in  the  front  and  rear  spar  web.  The 
shear  forces  applied  to  the  front  and  rear  spar  are: 
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Vh 


FS 


Fps  hFS2  +  hRS2 


and  Frs  — 


Vh 


RS 


hFS2  +  hRS2 


(D.2) 


Where  V  is  the  shear  force  applied  at  the  section  and  hFSand  hRS  are  the  height  of  the  front  and 
rear  spars  respectively. 


The  shear  flow  in  the  web  due  the  shear  force  is: 

Fpc  Fnc 

Qv,fs  =  7  and  Qv,rs  =  7  (D.3) 

nFS  nRS 

The  torque  creates  also  a  shear  flow  around  the  wing  box  in  addition  to  the  above  shear  flow  in 
the  spar  web. 

The  shear  flow  due  to  the  torque  is: 


T 


(D.4) 


Where  T  is  the  worst  value  of  the  torque  at  the  section  and  A  is  the  area  of  the  wing  box  cross- 
section. 

Because  of  the  orientation  of  the  shear  flow  in  the  spar  webs,  different  values  of  the  torque  have 
been  used.  As  presented  in  the  Figure  35,  the  maximum  positive  torque  value  is  used  for  the  front 
spar  whereas  the  maximum  negative  torque  is  applied  to  the  rear  spar. 

^  V  +ve 

\ 

Front  spar 

\ 


Figure  32:  Simplified  shear  flow  diagram  in  the  wing  box 

The  total  shear  flow  applied  to  the  front  and  rear  spar  webs  are: 

QfS  =  Qv,FS  +  Qt, positive  am ^  QrS  =  Qv,RS  ~  Qt, negative  (D.5) 


Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


63 


2013 


EOARD/AFRLl 


Cranfield 

I  UNIVERSITY 


Then,  the  different  sections  of  the  outboard  wing  were  sized  by  applying  the  previous  shear  flows 
to  the  laminate  in  CoALA  [29]  and  checking  the  failures  indices  and  strains. 

The  thicknesses  of  the  front  and  rear  spar  webs  are  presented  in  the  Table  6. 

Table  4:  Front  and  rear  spar  web  thicknesses  of  the  outboard  wing 


Front  Spar 

Section 

Web  Thickness  (mm) 

Rear  Spar 

Web  Thickness  (mm) 

5 

6.03 

2.88 

6 

5.50 

2.88 

7 

4.98 

2.62 

8 

4.45 

2.36 

9 

3.41 

1.83 

10 

2.36 

1.31 

11 

1.57 

1.05 

Regarding  the  inboard  sections,  the  increasing  torque  might  be  balanced  by  important  depth  of  the 
wing  box  and  the  drop  of  the  shear  force.  Since  the  front  spar  carries  most  of  the  shear  in  the 
outboard  wing,  the  thickness  of  the  web  has  been  kept  to  3.93  mm  for  the  sections  3  and  4  which 
are  the  most  critical  of  the  inboard  wing  in  terms  of  shear.  The  remaining  part  of  the  front  spar 
and  the  other  spars  of  the  inboard  wing  have  a  web  thickness  of  1.83  mm. 

The  spar  caps  are  the  components  making  the  link  between  the  spar  web  and  the  skin  surfaces. 
They  can  be  compared  to  stringers  as  they  support  mainly  bending.  The  layup  [+45/02/-45/90]s, 
same  as  the  skin/stringer,  and  T-section  beams  were  chosen  for  the  design.  Due  to  time 
constraints,  it  has  been  assumed  that  the  dimensions  are  double  the  stringers  dimensions. 
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D.2.4  Frames 


Heavy  Frames 

The  heavy  frames  were  sized  to  support  the  chordwise  bending.  They  support  compression  from 
the  spanwise  torque,  which,  in  the  chord  direction,  is  equivalent  to  the  bending.  The  caps  are 
designed  using  Euler  buckling  analysis.  The  layup  used  is  the  same  as  for  the  stringers,  [+45/02/- 
45/90]  s,  to  resist  the  compressive  loads. 

The  method  followed  is  iterative;  several  designs  have  been  tested  to  find  the  appropriate  design. 
First  of  all,  the  applied  loads  are  defined.  The  compressive  load  is  derived  from  the  spanwise 
torque  from  the  loading  actions: 


p  — _ v  i 

r  compressive  u 


(D.6) 


Where  T  is  the  ultimate  spanwise  torque  at  the  section,  A  the  total  wing  box  area  and  L  is  the 
length  of  the  frame. 

The  applied  stress  on  the  section  is  then  calculated: 

p 

_  _  r  compressive  /nk  ^ 

G applied  ~  1  (D.7) 

Aframe 


The  Euler  buckling  allowable  stress  is  computed  following  the  equation  [30]: 

_  n2EINA 

°cr~~mT 


(D.8) 


Where  E  is  the  tangent  modulus  of  the  material,  here  assumed  to  be  equal  to  the  elastic  constant 
E1;l  derived  from  CoALA,  INA  is  the  second  moment  of  inertia  about  the  neutral  axis,  A  is  the  area 
of  the  cross-section  of  the  beam  and  L  is  the  length  of  the  beam. 


Light  Frames 

The  light  frames  can  be  considered  as  “chordwise  stringers”.  They  are  designed  in  the  same  way. 
The  compressive  loads  applied  between  the  heavy  frames  are  calculated  the  same  manner  as  for 
the  heavy  frames.  The  layup  of  the  I-section  light  frames  is  also  the  same  as  the  stringers  and  the 
heavy  frames:  [+45/02/-45/90]s.  The  skin/light  frames  are  sized  with  the  automated  code  used  for 
the  skin/stringers  buckling  design.  The  thickest  skin  between  the  light  frame  and  stringer  skin 
design  is  finally  retained  for  the  final  section  skin  thickness. 
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Appendix  E.  Static  Finite  Element  Analysis 

FEA  (Finite  Element  Analysis)  is  a  computational  method  to  calculate  the  structural  behaviour 
under  specified  conditions.  The  FEA  is  commonly  used  in  research  departments  to  evaluate  the 
structure  strength  before  testing  it  with  expensive  experiments.  FEA  is  nowadays  a  key  step  in  the 
structural  design. 

In  the  framework  of  the  thesis,  the  wing  was  modelled  and  constrained  under  loads  with  the  pre¬ 
processor  Patran.  The  analyses  are  run  by  the  finite  element  calculator  Nastran. 

In  this  chapter,  the  static  linear  analysis  of  the  half  of  the  wing  is  presented.  Further  analyses  of 
the  structure  modes  and  dynamic  response  are  described  in  the  later  chapters. 

Because  FEA  can  be  a  time-consuming  task,  the  author  built  the  model  step  by  step  to  avoid 
numerous  problems.  A  mesh  sensitivity  analysis  has  been  conducted  to  evaluate  the  reliability  of 
the  results  of  the  static  analysis  in  function  of  the  mesh  size. 

E.l  CATIA  Surface  Model 

The  geometry  of  the  wing  was  constructed  in  the  software  CATIA.  The  skins,  spar  and  heavy 
frame  webs  and  ribs  were  represented  by  surfaces.  The  beams,  such  as  the  stringers,  spar  caps  and 
frame  beams  were  symbolised  by  lines.  All  the  surfaces  were  cut  by  the  spar,  rib,  frame  and 
stringer  planes  to  permit  an  easier  mesh  in  Patran.  The  surface  model  was  imported  into  Patran 
using  the  Parasolid  import  function  for  the  STEP  files  from  CATIA.  The  inboard,  outboard  and 
device  part  were  separated  into  different  groups  to  allow  later  changes  in  the  device  design  for 
example.  The  final  surface  model  counts  1922  surfaces  for  half  of  the  aircraft. 
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Figure  33:  Surface  model  from  CATIA 

E.2  Mesh 

A  mesh  sensitivity  analysis  has  been  conducted  to  evaluate  the  influence  of  mesh  element  size  on 
the  results.  Three  different  meshes  were  built  in  Patran,  by  increasing  the  number  of  mesh 
elements  by  approximately  two  between  each  mesh.  The  number  of  elements  and  nodes  is  given 
in  the  Table  7. 


Table  5:  Mesh  nodes  and  element  numbers 


Mesh  1 

Mesh  2 

Mesh  3 

Nominal  element  edge  length  (mm) 

200 

150 

100 

Number  of  Elements 

19169 

29657 

63130 

Number  of  Nodes 

11010 

19334 

47106 

The  meshes  of  the  surfaces  use  mainly  Quad4  elements  with  the  IsoMesh  function.  In  the  other 
cases,  Tria3  elements  were  used  and  paver  and  hybrid  options  were  chosen  when  the  surfaces 
were  not  rectangular  and  impossible  to  be  meshed  with  IsoMesh  function. 

To  ensure  that  the  elements  are  all  connected  together,  mesh  seeds  were  used  for  the  edges. 
Equivalence  of  the  nodes  has  been  completed  and  free  edges  have  been  checked.  The  verification 
function  of  the  elements  in  Patran  was  used  to  check  the  dimensions  of  the  Quad  and  Tria 
elements  (aspect  ratio<5,  skew  angle>30  for  quad,  skew  angle>10,  taper<0.5  and  warp  angle<0.05 
for  quad). 
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The  stringers,  frames  and  spar  caps  were  meshed  by  bar  elements  along  the  lines  of  the  surface 
model. 


Figure  34:  Mesh  1 


E.3  Properties 

The  material  properties  were  input  into  the  software.  In  the  first  approach,  an  equivalent 
aluminium  alloy  was  used  to  check  that  the  meshes  have  not  got  any  errors.  Then  the  properties  of 
the  8552/IM7  UD  composite  prepreg  were  entered  into  the  software  as  a  2D  orthotropic  material. 
The  layups  of  the  laminates  were  defined  for  the  components,  following  the  choices  done  during 
the  initial  sizing.  The  shell  properties  for  the  skins,  ribs,  frame  webs  and  spar  webs  as  well  as  the 
beams  dimensions  for  the  stringers,  frames  and  spar  caps  were  attributed  to  the  mesh  elements. 

E.4  Boundary  Conditions  and  Loads 

The  model  input  in  Patran  is  the  half  of  the  aircraft.  To  represent  the  boundary  conditions  at  the 
root  section,  the  nodes  on  the  edges  at  the  centreline  were  fixed,  constraint  in  both  translation  and 
rotation. 

The  loads  applied  on  the  structure  were  derived  from  the  previous  aerodynamic  results.  Unlike  the 
initial  sizing  where  shear,  bending  and  torque  need  to  be  calculated  to  size  the  components,  the  FE 
model  only  needs  the  aerodynamic  loads.  Indeed,  the  results  of  loads  into  bending  moments  are 
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part  of  the  calculation  of  the  software.  The  mass  of  the  structure  is  also  taken  into  account  by  the 
application  of  an  inertial  load. 

The  worst  load  case  in  terms  of  bending  was  used  for  the  analysis.  This  is  the  case  at  sea  level  at 
empty  weight  for  a  static  gust  load  factor  of  3.81  g  (gust  case  2).  The  aerodynamic  loads  applied  in 
the  model  are  derived  from  the  lift  distribution  obtained  during  the  aerodynamic  analysis.  The  lift 
and  the  pitching  moment  are  distributed  along  the  span  by  idealising  it  in  punctual  loads.  The 
loads  are  applied  on  the  nodes  located  at  the  aerodynamic  centre  of  the  frame  and  rib  webs. 

The  inertia  of  the  components  is  considered  by  the  application  of  an  inertial  load  of  3.81  g  to  the 
whole  structure. 

Figure  38  presents  the  total  load  distribution  and  boundary  condition  applied  on  the  model. 


Inertia 


Figure  35:  Loads  and  boundary  conditions  applied  on  the  model 

E.5  Static  Analysis  Results 

Once  the  model  was  completely  built,  a  static  linear  analysis  was  performed  by  the  solver  Nastran 
for  each  of  the  three  meshes. 

E.5.1  Results  for  the  Mesh  1 

Figure  39  presents  the  displacement  of  the  wing  under  limit  loads.  The  wing  tip  reaches  a 
maximum  displacement  of  2.56  m.  This  value  demonstrates  the  flexible  characteristic  of  the  flying 
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wing,  but  it  needs  to  be  put  in  perspective  with  the  fact  that  it  is  obtained  for  the  load  factor  of 


3.81g. 
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Figure  36:  Displacement  of  the  structure  under  limit  loads 

The  failure  indices  (FI)  and  the  strains  of  the  members  preliminary  sized  were  computed  to  verify 
that  they  are  under  the  design  limits.  The  Twai-Wu  failure  criterion  was  used  for  the  calculation  of 
the  FI. 

It  can  be  seen  in  Figure  40  that  the  FI  do  not  exceed  the  value  of  0.56  in  the  upper  skin.  The 
strains  are  plotted  as  minimum  principal  strain  because  the  upper  skin  supports  mainly 
compression  loads.  The  maximum  magnitude  of  the  strains  reaches  3320  ps  which  is  under  the 
limit  of  3500  ps  established  at  the  beginning  of  the  study  (Figure  41).  It  can  be  observed  that  the 
main  strains  and  stresses  are  located  around  the  kink,  in  the  outboard  part,  where  the  loads  are 
high  and  the  structure  thin. 
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Figure  37:  Upper  skin  FI 
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Figure  38:  Upper  skin  strains 

The  lower  skin  is  under  tensile  loads.  The  FI,  illustrated  in  the  Figure  54,  stay  under  one  with  a 
maximum  value  of  0.81.  However  the  strains  are  critical  for  the  design  of  composite  structure. 
The  design  limit  of  3600  ps  is  exceeded  at  the  kink  with  a  maximum  strain  of  5200  ps  (Figure 
43).  The  stresses  at  this  concentration  area  must  be  reduced  to  avoid  too  high  strains  for  damage 
tolerance  considerations.  After  verification  that  this  concentration  of  loads  is  not  due  to  the  mesh, 
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the  thicknesses  of  the  components  in  this  region  were  increased  to  have  the  strains  under  the 


design  limit. 
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Figure  39:  Lower  skin  FI 
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Figure  40:  Lower  skin  strains 

The  spars  FI  and  strains  are  also  checked.  Excluding  the  concentration  points  in  the  inboard  spars 
due  to  the  application  of  punctual  loads  in  the  model,  the  same  issue  as  the  lower  skin  is 
encountered  for  the  rear  spar  at  the  kink.  Although  the  maximum  FI  is  0.86  at  this  location  (Figure 
44),  the  strains  reach  the  upper  value  of  5200  ps.  In  the  same  approach  as  for  the  lower  skin,  the 
thickness  of  the  rear  spar  was  increased. 
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Figure  42:  Outboard  spar  strains 


E.5.2  Mesh  Sensitivity  Study 

This  is  a  primary  concern  to  ensure  that  the  results  obtained  by  the  FE  model  are  reliable.  The 
mesh  sensitivity  analysis  is  used  to  verify  that  the  results  are  not  a  function  of  the  spatial 
discretisation  of  the  model.  The  results  presented  above  from  the  first  mesh  can  be  compared  to 
the  results  obtained  with  the  two  finer  meshes.  The  deflection  of  the  wing,  the  FI  and  strains  for 
the  components  were  calculated  for  each  of  the  meshes.  The  results  are  presented  in  the  Table  8. 
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Table  6:  Comparison  of  the  results  between  the  meshes  (worst  values) 


Mesh  1 

(coarse) 

Mesh  2 

Mesh  3 

(finest) 

Deflection  (limit  loads) 

2.56  m 

2.56  m 

2.55  m 

Upper  skin  FI 

0.51 

0.52 

0.53 

Upper  skin  strains 

3320  pe 

3160  pe 

3520  pe 

Lower  skin  FI 

0.81 

0.79 

0.91 

Lower  skin  strains 

5200  pe 

5050  pe 

4700  pe 

Outboard  rear  spar  FI 

0.86 

0.86 

0.92 

Outboard  rear  spar  strains 

5200  pe 

5050  pe 

5250  pe 

In  the  Figure  46,  it  can  be  seen  that  the  results  from  the  three  meshes  are  quite  convergent.  If  mesh 
3,  the  finer  mesh,  is  taken  as  the  reference,  the  maximum  difference  between  the  values  is  13%  for 
the  lower  skin  FI  between  the  meshes  2  and  3. 

Therefore,  the  results  obtained  with  the  mesh  1  can  be  considered  reliable  enough  considering  the 
assumptions  done  in  terms  of  modeling  of  the  loads.  This  mesh  will  be  used  for  the  further  finite 
element  analyses  which  demand  a  higher  amount  of  calculations  and  computational  resources. 
The  optimization  process  also  limits  the  number  of  mesh  elements. 


Figure  43:  FI  results  in  function  of  the  number  of  mesh  elements 
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E.5.3  Update  of  the  Structural  Component  Dimensions 

As  presented  earlier,  the  strains  in  the  lower  skin  and  outboard  rear  spar  at  the  kink  are  exceeded. 
In  order  to  reduce  it  to  the  design  limit,  the  thicknesses  of  the  components  of  the  design  section  5 
have  been  increased  progressively.  The  dimensions  of  the  lower  skin  and  rear  spar  are  given  in  the 
Table  9. 


Table  7:  Updated  dimensions  of  the  lower  skin  and  rear  spar 


Lower  skin 

Rear  spar 

Previous  thickness  (mm) 

4.45 

2.89 

Updated  thickness  (mm) 

6.29 

5.50 

With  these  new  dimensions,  the  highest  strains  are  dropped  from  5200  ps  to  3570  ps  as  shown  in 
Fig.47,  which  meets  the  specified  design  requirement.  The  wing  tip  deflection  of  this  updated 
design  reaches  2.4  m  under  limit  load. 

Palran  2010.1 .2  64-Bit(MD  Enabled)  10Aug-l2  1 1  10:41  3.66-003 

Fringe  Uwtjoads.  A1 7  Static  Subcase.  Strain  Tensor. .  Max  Pnnopal.  68  ol  68  layers  ( Maximum)  333-003 

Deform  Lfnrt_loads.  A1 7:Sta1ic  Subcase.  Displacements.  Translational.  3  09^3 

286-003 


Max 

Figure  44:  Strains  in  the  lower  skin  and  outboard  rear  spar  (updated  dimensions) 
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Appendix  F.  Additional  Gust  Analysis  and  BDF  Code  in  NASTRAN 

F.l  Gust  for  sweepback  shaft  parallel  to  front  spar 

In  this  case,  the  shaft  was  set  to  parallel  to  the  front  spar  for  checking  the  difference  in  gust 

alleviation.  It  was  easy  to  find  that  the  efficiency  of  sweepback  shaft  configuration  is  not  as  high 
as  the  straight  shaft. 


Table  8  gust  response  with  sweepback  shaft  PGAD 


Case 

Spring 

stiffness 

(Nm/rad) 

Wing  tip 
Disp. 

(m) 

Disp. 

Reduction 

Bending 

Moment 

(KNm) 

BM 

Reduction 

PGAD 

Relative 

Twist 

angle(°) 

Initial 

design 

/ 

3.07 

/ 

5600 

/ 

/ 

a=-0.7 

5.8E4 

2.81 

8.5% 

5195 

7.2% 

6.5 

a=-0.5 

5.8E4 

2.97 

3.3% 

5446 

2.8% 

4.2 

4000 

3000 

2000 

1000 

0 

-1000 

-2000 


A 

a=-0.7,sweeo|)ack  shaft 
h  a=4jf.5,sweepback  shaft 

- without  PGAD 

; 

Time,  s 

0 

Vr  4  ' 

Figure  48  wing  tip  displacement  response 


Figure  49  PGAD  relative  twist  angle 


Figure  50  wing  bending  moment  response 
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Figure  48-50  presented  the  analysis  results  in  terms  of  wing  tip  displacement,  PGAD  relative 
twist  angle  and  wing  root  bending  moment.  The  device  twist  angle  reduced  dramatically 
comparing  with  the  straight  shaft  so  that  there  was  just  a  slightly  alleviation  about  2.8%  for  wing 
root  bending  moment. 

The  following  codes  are  the  core  part  of  Nastran  input  file.  Wing  structure  with  PGAD  were 
introduced  in  previous  appendices  in  details  and  not  included  here. 


F.2  BDF  Code  used  in  NASTRAN 


$■ 


$ 


Control  code- 


$  Dynamic  Gust  Analysis 
SOL  146 
TIME  600 
CEND 

$  Direct  Text  Input  for  Global  Case  Control  Data 
TITLE  =  MSC.Nastran  Aeroelastic  job  created  on  06-Dec-12  at  14:24:34 
ECHO  =  NONE 
MAXLINES  =  999999 
AECONFIG  =  aeroPGAD 
SUBCASE  1 
$  Subcase  name  :  test 
SUBTITLE=Default 
METHOD  =  1 

SDAMP  =  2000  $  STRUCTURAL  DAMPING  (2  PERCENT) 

GUST  =  1000  $  AERODYNAMIC  LOADING  (1  -cos  GUST) 

DLO AD  =  1001  $  REQUIRED 

FREQ  =  40  $  FREQUENCY  LIST 

TSTEP  =  41  $  SOLUTION  TIME  STEPS  (1  PERIOD) 

SPC  =  2 

AESYMXZ  =  Symmetric 
AESYMXY  =  Asymmetric 

$ - GLOBAL  DEFINATION - $ 

$  Direct  Text  Input  for  Bulk  Data 
PARAM  POST  0 
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PARAM 

WTMASS  1. 

PARAM 

SNORM  20. 

PARAM 

PRTMAXIM  YES 

PARAM 

GUSTAERO 

-1 

PARAM 

MACH  0.3 

PARAM 

Q  4.965-3 

PARAM 

LMODES  11 

EIGRL 

1 

10 

$ - GUST  information - $ 

$  ID  TYPE 

TABDMP1  2000  +TABDMP 

$  FI  G1  F2  G2  "ENDT" 

+TABDMP  0.  .03  10.  .03  ENDT 

$  SID  DLOAD  WG  X0  V 
GUST  1000  1001  0.189  -25000.  101000.0 

$  SID  DAREA  DELAY  TYPE  TID 
TLOAD1  1001  1002  1003 


$  DAREA  DEFINES  THE  DOF  WHERE  THE  LOAD  IS  APPLIED  AND  A  SCALE  FACTOR. 

$ 

$  SID  P  C  A 

DAREA  1002  90000  3  0. 

TABLED  1  1003  +TAB1 

$  XI  Y1  X2  Y2  X3  Y3  X4  Y4 

+TAB1  0.27692  0.02446  0.35385  0.09548  0.43076  0.20608  0.50769  0.34547  +TAB2 
+TAB2  0.58462  0.49994  0.66154  0.65447  0.73846  0.79386  0.81539  0.90446  +TAB3 
+TAB3  0.89231  0.97548  0.96923  1.00000  1.04615  0.97548  1.12308  0.90446  +TAB4 
+TAB4  1.2  0.79381  1.27692  0.65447  1.35385  0.49994  1.43077  0.34547  +TAB5 

+TAB5  1.50769  0.20608  1.58462  0.09548  1.66154  0.02446  1.73846  0  +TAB6 


+TAB6 

ENDT 

$  SID  N 

DT  NO 

FREQ1 

40 

0. 

.125 

TSTEP 

41 

600 

.01 

$ - spring  definition - $ 

$  Elements  and  Element  Properties  for  region  :  rigidspring_ux 
PELAS  108  1.+10 

$  Pset:  "rigidspring_ux"  will  be  imported  as:  "pelas.108" 

CELAS1  60000  108  90000  1  90001  1 
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$  Elements  and  Element  Properties  for  region  :  rigidspring_uy 
PELAS  109  1.+10 

$  Pset:  "rigidspring_uyn  will  be  imported  as:  "pelas.109" 

CELAS1  60001  109  90000  2  90001  2 

$  Elements  and  Element  Properties  for  region  :  rigidspring_uz 
PELAS  110  1.+10 

$  Pset:  "rigidspring_uz"  will  be  imported  as:  "pelas.110" 

CELAS1  60002  110  90000  3  90001  3 

$  Elements  and  Element  Properties  for  region  :  rigidspring_Rx 
PELAS  111  1.+10 

$  Pset:  "rigidspring_Rx"  will  be  imported  as:  "pelas.lll" 

CELAS1  60003  111  90000  4  90001  4 

$  Elements  and  Element  Properties  for  region  :  rigidspring_Ry 
PELAS  112  1.+7 

$  Pset:  "rigidspring_Ry"  will  be  imported  as:  "pelas.112" 

CELAS1  60004  112  90000  5  90001  5 

$  Elements  and  Element  Properties  for  region  :  rigidspring_Rz 
PELAS  113  1.+10 

$  Pset:  "rigidspring_Rz"  will  be  imported  as:  "pelas.113" 

CELAS1  60005  113  90000  6  90001  6 

$ - AERODYNAMICS - $ 

$  MKAER02  card 


$  Mach-Frequency  Pair  .MRG_MK_2 


MKAER02  .265 

.2 

.265 

.4 

.265 

.6 

.265 

.8 

MKAER02  .265 

1. 

.265 

1.2 

.265 

1.4 

.265 

1.6 

MKAER02  .265 

1.8 

.265 

2. 

.265 

2.2 

.265 

2.4 

MKAER02  .265 

2.6 

.265 

2.8 

.265 

3. 

.265 

3.2 

MKAER02  .265 

3.4 

.265 

3.6 

$ 

$  Aeroelastic  Model  Parameters 
PARAM  AUNITS  1. 

$ 

$  Global  Data  for  Steady  Aerodynamics 

$ 

$  A  half-span  model  is  defined 

$ 

AERO  0  90000.  6260.  1.227-12 

AEROS  0  0  6260.  62280.  1.949+08 

$ 
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$  Flat  Aero  Surface:  Device 


Cranfield 
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PAEROl  102001 

CAEROl  102001  102001  0  2  4  1 

10096.  29140.4  0.  4563.  10920.  31140.  0.  4563. 

$  Flat  Aero  Surface:  outboard_wing 


PAEROl  101001 

CAEROl  101001  101001  0  15  4  1 

-768.  10795.  0.  4563.  10096.  29140.  0.  4563. 


$  Flat  Aero  Surface:  inboard_wing 


PAEROl  100001 

CAEROl  100001  100001  0  5  4  1 

-7000.  0.  0.  14706.  -768.  10795.  0.  4563. 

$ 

$  Surface  Spline:  spline_rib_nodes_l 

$ 

SPLINE4  1  100001  1  1  IPS  BOTH 

10  10 

AELIST  1  100001  100002  100003  100004  100005  100006  100007 

100008  100009  100010  100011  100012  100013  100014  100015 

100016  100017  100018  100019  100020  101001  101002  101003 

101004  101005  101006  101007  101008  101009  101010  101011 

101012  101013  101014  101015  101016  101017  101018  101019 

101020  101021  101022  101023  101024  101025  101026  101027 

101028  101029  101030  101031  101032  101033  101034  101035 

101036  101037  101038  101039  101040  101041  101042  101043 

101044  101045  101046  101047  101048  101049  101050  101051 

101052  101053  101054  101055  101056  101057  101058  101059 

101060 


SET1  1 

1 

4 

19  43  79 

123 

176 

236 

306 

385 

473  570  676  684  791 

808 

913 

928 

1033 

1048 

1153 

1168  1273 

1288 

1393 

1408 

1513 

1528 

1633 

1648 

1753 

1768 

1873 

1888 

1993 

2008 

2113 

2128 

2233 

2248 

2352 

2366 

2477 

2579 

2672 

2757 

2833 

2899 

2956 

3004 

3044 

3078 

3099 

17021 

17162 

17603  17872  18594  19005  19012  19504  20698  21386 
SPLINE4  2  102001  2  2  IPS  BOTH 
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10  10 

AELIST  2  102001  102002  102003  102004  102005  102006  102007 

102008 

SET1  2  40537  40561  40585  40743  40755  40767 
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